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ABSTRACT 


Nisheeth Patel, Doctor of Philosophy, 1983 

Major: Engineering, Department of Aerospace Engineering 

Title of Dissertation: A Fully Vectori 2 ed Numerical Solution of the 

Incompressible Navier-Stokes Equations 
Directed by: Joe F. Thompson 

Pages in Disseration: 160 Words in Abstract: 255 

Abstract 

A. vectorizable algorithm is presented for the implicit finite 
difference solution of the incompressible Navier-Stokes equations 
in general curvilinear coordinates. The unsteady Reynolds averaged 
Navier-Stokes equations solved are in two-dimension and non-conservative 
primitive variable form. A two-layer algebraic eddy viscosity turbu- 
lence model is used to incorporate the effects of turbulence. Two 
momentum equations and a Poisson pressure equation, which is obtained by 
taking the divergence of the momentum equations and satisfying the contin- 
uity equation, are solved simultaneously at each time step. An elliptic 
grid generation approach is used to generate a boundary-conforming 
coordinate system about an airfoil. The governing equations are express- 
ed in terms of the curvilinear coordinates and are solved on a uniform 
rectangular computational domain.. A checkerboard SOR, which can effect- 
ively utilize the computer architectural concept of vector processing, is 
used for iterative solution of the governing equations. The method 
is applied to the cases of an 18% thick NACA 66^018 airfoil at zero 
degree angle of attack for chord Reynolds number range of 1000-40,000. 

v.i 




The effects of various boundary-conforming coordinate systems, arti- 
ficial viscosities, smoothers, down-stream boundary conditions, initial 
guesses and number of iterations during the acceleration phase on the 
solution of the flow field are studied. Numerical results are given in 
terms of surface pressure distributions and velocity vector fields at 
selected times. Computed steady-state results are compared with 
experimental data. On the CDC CYBER-203 computer, the algorithm 
demonstrated a factor of about ll improvement over a CDC-7600 scalar 
version of the code. 
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Chapter I 
INTRODUCTION 

Computational Fluid Dynamics (CFD) has made a significant contri- 
bution in the recent development of aerospace vehicles. Practical 
aerodynamics, which controls the design of flight vehicles, is 
essentially about complex flow at high Reynolds number past arbitrary 
configurations. The governing equations, which describe physical 
features of such a flow, are non-linear partial differential 
equations - the Navier-Stokes equations. Simplification of these 
governing equations will limit the application. In the past, experi- 
mental fluid dynamics lias played an important role, however, with 
the breakthrough in solving non-linear partial differential equations 
and high speed computation, CFD has risen to complement the role of 
experimental fluid dynamics. 

Perhaps the foundation of modern fluid dynamics was laid by Prandtl 
when he first presented boundary layer theory in 1904. However, it was 
not recognised until 1920 when Prandtl presented insight on separation. 

A review of classical fluid dynamics has been presented by Goldstein 
[1]. The foundation stone for CFD was probably laid by Courant, 
Friedricks and Levy [2] with the introduction of the numerical stability 
condition for the solution of hyperbolic equations, known as CFL 
condition. In the 1960's the finite difference methods, such as 
Lax -Wendrof'f [3] type and MacCormack's [4] explicit methods, were 
developed for solving the Euler equations in conservation law form. 

Since then, with rapid increase in computer speed and computer memory, 

CFD has developed sufficiently to become established as a discipline. 
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Several recent surveys on CFD development and future of CFD have been 
presented in references [5], [6], [7], and [8], 

Although in most cases CFD does offer the potential of obtaining 
complete information about complex flows without experimentation, 
it has its own difficulties. The essential areas to be considered 
before solving the Navier-Stokes equations are grid generation, 
algorithm, turbulence model, and computer. Significant improvement 
in any of these areas will enhance computational efficiency and 
accuracy of the solution. 

Many cases of practical interest contain an arbitrary domain. 

Since the boundary is not aligned with the grid when using a cartesian 
coordinate for an arbitrary region, it requires the use of inter- 
polation formulas near the boundary. The Imposition of boundary 
conditions with a complicated computational region having irregular 
boundaries is a primary difficulty with the cartesian coordinate 
system. Moreover, the Navier-Stokes equations and their boundary 
conditions are such that the viscous effects are confined to a very 
thin region immediately adjacent to the solid boundaries. Although 
the region is quite thin it produces considerable effects on the total 
solution of the flow field. In addition, the stability conditions, 
iterative convergence and truncation errors of the numerical algorithm 
employed may be adversely affected [9]. Since a cartesian grid has 
limited applications, the recent trend has been to use a boundary- 
conforming coordinate system. The boundary-conforming coordinate 
systems are defined as those which possess constant coordinate lines 
coincident with all boundaries of the physical plane, which in turn 
correspond to a rectangular grid with square cells in the transformed 
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plane. The governing continuum equations are derived on the rectang- 
ular grid in the transformed plane. An elliptic grid generating 
system developed by Thompson, et. al. , [101 is capable of generating 
a boundary-conforming coordinate system. The main attraction of this 
approach is flexibility, automation and a moderate degree of control 
by the user. The recent surveys on grid generation techniques and 
applications are presented by Thompson in references [11], [12], and 

[13]. 

In the past decade, numerical algorithms used in simulation of 
fluid flows have improved substantially. Explicit algorithms are 
simple. However, restriction on the time step imposed by stability 
considerations is a main disadvantage of these schemes. Increased 
interest in implicit schemes led to the development and use of efficient 
algorithms such as those due to Briley and McDonald [14], Beam and 
Warming [15], MacCormack's rapid solver [16], and the hybrid MacCor- 
mack's scheme [17], The mathematical reviews of these developments 
are presented by Lomax in references [18] and [7]. 

Practical computations involve numerical simulation of turbulence 
to provide more understanding of physical phenomena. Several basic 
algebraic, one-equation and two-equation models have been developed 
and used to analyze turbulent flows. A comprehensive review on 
turbulence modelling has been presented by Marvin in reference [19]. 

In the past decade computer speeds and computer memories have 
increased at a significant rate. The development of supercomputers, 
with memory measured in million words and calculation rate in mflops 
(million floating point operations per second) , such as the ILLIAC-IV 
(16 million words, 25 mflops), CYBER-203 (1 million words, 20 mflops), 
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CRAY-1 (1 million word, 30 mflops) , CYBER-205 (4 million words, 80 mflops) 
and C.RAY-1S (4 million words, 30 mflops) have dramatically increased 
the capabilities of CFD and reduced the cost of computation. For 
example, the CRAY-1S can perform 100,000 calculations for less than a 
penny. The. CRAY -2, which is still in the design stage, will approximately 
double the speed of the CRAY-1S with possible further reduction in 
calculation cost. The supercomputers are used in areas such as 
aeronautical engineering, nuclear research, weather forecasting and 
other military and civilian applications. The CYBER-200 and CRAY-1 
series computers are vector processors and use pipe-line architecture 
to increase, the calculation rate. Levine presented an introduction 
to supercomputers and their architecture in reference [20], and tech- 
nical information can be found in reference [21]. A computational 
algorithm developed with the supercomputer architecture in mind can 
effectively use its computational capabilities and hence reduce the 
run costs significantly. 

The present study will be focused in the areas of algorithm 
development and the use of a supercomputer for the numerical solution 
of incompressible Navier-Stokes equations. The areas of grid genera- 
tion and turbulence, modelling will be addressed as essential elements 
required for simulation of the flow. The information about grid 
generation techniques and the turbulence model used in this study 
can be found in appropriate references. However they will be presented 
in detail for completeness. 
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The unsteady and steady incompressible Navier-Stok.es equations 
have been successfully used by many researchers to simulate the flow 
field of different characteristics. The basic formulations used in 
most work are velocity-vorticity, streamfunction-vorticity , stream- 
function-biharmonics and primitive variables. Cebeci has reviewed 
the last three formulations in reference {22]. 

The viscous incompressible flow past an airfoil has been subjected 
to several numerical attacks in the past decade. The pioneering 
computations of laminar, incompressible two-dimensional flows about 
an airfoil are summarized in [23] and [24]. Hodge [25] used the 
optimized boundary conforming coordinate system for the laminar flow. 

In recent years, turbulent flows have been of increasing interest 
to researchers. Sugavanam [26] obtained the solution for flow 
past a Joukowski airfoil using velocity-vorticity formulation. Hegna 
[27] used primite variables for flow past a NACA 0012 airfoil. 

Bernard [28] employed the Approximate Factorization technique for 
NACA-66 a -0018 airfoil section. Moif.ra [29] simulated three dimensional 
turbulent flow around an airfoil. Lately, Freeman [30] used an adaptive 
grid approach for dynamic coupling of the grid and flow field solution. 
Experimental investigations have been reported by Mueller in reference 
[31]. 

1 . 2 Research Objectives and Outline 

The objective of this effort is to develop a vectorized computer 
code for viscous turbulent, two-dimensional incompressible flow past 
an airfoil using an implicit finite-differencing scheme (Backward- 
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Time, Central-Space) . On the pipeline computers, such as the CYBER- 
200 series, it is desirable to work with very long vectors for efficient 
use of its vector processing capabilities. Explicit methods are sim- 
ple and can be easily vectorized since the entire grid can be considered 
as a long vector,. However, the major disadvantage of explicit schemes 
lies in the time step restrictions imposed by stability considerations. 
Implicit schemes are frequently unconditionally stable and usually 
employ an iterative method such as SOR (Successive Over Relaxation) . 

Since application of SOR on a vector machine results in either the 
inefficient use of its vector processing abilities or the necessity 
to shift to slower scalar operations, the checkerboard SOR algorithm 
(Chapter VII) will be used for the iterative solution of the governing 
equations on a vector processor. Also the effects of the algorithm, 
smoothers, grid, various forms of articifial viscosity and some boundary 
conditions on the solution will be investigated for the specific 
problem under study. The comparisons will be made with the available 
experimental results in Chapter VIII. 

The present research effort is carried out in the following 
frame-work. The governing equations are the two-dimensional, 
incompressible Reynolds-averaged Navier-Stokes equations, written in 
non-conservative form in terms of primitive variables. A Poisson 
equation for the pressure is obtained by taking the divergence of the 
momentum equations. The two -layer algebraic turbulence model of 
Baldwin and Lomax [32] is used to calculate the eddy viscosity for 
the Reynolds-averaged equations. Dirichlet boundary conditions are 
imposed on the freestream boundary. On the downstream boundary , extra- 
polation boundary conditions on the velocity and Dirichlet boundary 
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condition on the pressure are imposed. The boundary conditions imposed 
on the airfoil surface are obtained employing no-slip conditions for 
the velocity and by setting the normal derivative of the pressure equal 
to zero on the boundaries. A linear gradual start is used to accelerate 
the flow from rest to its final freestream velocity. An implicit 
finite-differencing scheme, obtained using backward-time central-space 
approximations for the governing equations in the transformed plane is 
used to obtain flow field solutions. At each time step, the three govern- 
ing equations are solved simultaneously, for u, v, and p. 
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THE BOUNDARY-CONFORMING CURVILINEAR COORDINATE SYSTEM 
2 . 1 The Boundary-Conforming Concep t 

A coordinate system can have a significant influence on the 
numerical solution of hosted partial differential equations. For 
many cases of practical interest, the irregularities present in the 
boundary geometry will limit the use of the Cartesian coordinate 
system in finite difference flow field simulation. A cartesian coor- 
dinate system under such circumstances will require interpolation near 
the body boundary to implement the boundary conditions. In the boundary 
conforming coordinate system grid lines coincide with the body boundary 
thus yielding a degree of simplicity in the implementation of the bound- 
ary conditions. Also, the use of boundary conforming grids in the solu- 
tion of partial differential equations in domains surrounding arbitrary 
geometrical boundary shapes will give a well-ordered ssystem of algebraic 
difference equations compatible with the algorithms which can efficient- 
ly use the vector-processing computers. Various possible approaches 
such as conformal mapping, transfinite mapping, algebraic and elliptic 
equations have been successfully employed to generate body-conforming 
curvilinear coordinate systems. A comprehensive survey to these 
techniques and applications have been given by Thompson in references 
[11] and [13] and by Thompson, Wars! and Mastin in reference [12]. 

The feasible and systematic way of generating an appropriate 
body-conforming coordinate system should consider constraints that 
arc?, needed for a given problem. Preferably, the grid should be generated 
in an automatic manner with the elements of control within the mesh 
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generation process. Also the user should have an acceptable degree of 
control over grid smoothness, skewness and stretching. Hosted algo- 
rithms are usually sensitive to the grid smoothness, skewness and 
stretching and general reasons for this effect include the following: 

The coefficients of the transformed partial differential equations 
depend on the derivatives of the functions defining the. coordinate 
system thus smoothness of the grid will have considerable effects on 
the accuracy of the solution. The local truncation error increases 
with departure from orthogonality. Also, the use of algebraic turbulence 
models demand near-orthogonality at the boundary for consistent modelling 
of turbulent flows. For a fixed number of grid points, the clustering 
of grid points in the region of large gradient should reduce the error and 
improve solutions. Moreover, some algorithms require a grid generating 
procedure that can be dynamically coupled to the physical solution 
properties to enhance accuracy and efficiency of the numerical results. 

The grid generated using conformal mapping techniques have been 
used by several investigators. The main advantage is that it allows 
greater control by the user. The main disadvantage of this method is 
the lack of flexibility and automation. An elliptic system can generate a 
grid in an automatic manner with a moderate degree of control by the 
user.. It can be extended to three-dimensions and a .ptive grid 
procedures. 

2 . 2 Elliptic Systems 

An elliptic grid generating system proposed by Thompson, Thames 
and Mastin [10] and successfully used by many investigators is capable 
of generating a suitable grid for the present study. The elliptic 
grid generating system is less susceptible to grid overlapping and 
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can be subjected to a variety of grid control procedures to obtain 
desired grid characteristics as discussed in the previous section. 

Numerical grid generation usually involves transformation of 
the physical domain of interest into a geometrically simple compu- 
tational domain, such as a single rectangular domain. The solution 
of grid generation equations in the computational domain produces 
the corresponding grid in the physical domain. The physical space 
defined by Cartesian coordinates x and y is mapped onto the compu- 
tational space through the mapping functions 

£ = £(x,y) 
n = n(x,y) 

by making the inner, outer, lower downstream and upper downstream 

boundaries coincide with n=n. n = n , £ = £ . and £ = £ » 

min, max’ ^min max 

respectively. The extremum principle, ensures occurence, of extrema 
only on the boundaries and hence overlapping of grid lines can be 
avoided. 

The topological correspondence in a C-type grid about a 2-D air- 
foil may be better understood withthe help of Figure 1. The boundary 
n * r l m ^ n is mapped onto the inner boundary - T - r,j containing the 

branch cuts, and the airfoil £ = £ , and £ = £ correspond to the 

min max ^ 

downstream section r^ L and r 3u respectively, £ increasing clockwise 
around the airfoil. The n = constant: family of lines form open curves 
resembling the letter C. The q = q^^ boundary is mapped on to the 
outer freestream boundary T 2 - 

The elliptic grid generation method of Thompson et.al. [10] 
permits any desired distribution of £ and n on the boundaries. The 


( 2 . 1 ) 

( 2 . 2 ) 
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inherent smoothness of solutions of elliptic systems is well recognized, 
and hence they are less adaptable to propagation of boundary slope 
discontinuities into the field. The choice of the elliptic system is 
further reinforced by the ability of the inhomogenous terms in Poisson's 
equation to control coordinate line, spacing with respect to a curve or a 
point within the field. The chosen grid generating system has the 
following form: 

(2.3) 


a 


+ Eyy • 7 PU ' n) 


n + n = -t- Q(5»n) 

xx yy 2 ^ ’ 


(2.4) 


A desired form of the control functions P and Q makes it possible to 
concentrate lines in regions of the field. An interchange of dependent 
and independent variables enables one to perform all computation in the 
transformed field (Appendix A). The generating system in the transformed 
field becomes: 


ax -26x_ + yx “ -(ax P + yx Q) 

55 5n nn 5 n 


ay 55" 2f5y 5n + Yy nn = " (ay 5 P + Yy n Q) 


(2.5) 

( 2 . 6 ) 


The transformed equations are solved in the rectangular §n- plane and 
the Dirichlet boundary conditions are specified for x and y by the 
known shape of boundaries. The coefficients of equations (2.3) - (2.6) 
are functions of the transformation and are defined by 

(2.7) 


2 , 2 
a = x + y 
n n 


(5 = + y 5 y n 

2.2 
Y - x c + y {; 


( 2 . 8 ) 

(2.9) 
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J ■ Vn " Vc 


( 2 . 10 ) 


where J is the Jacobian of the transformation. 

A great deal of simplification in computation results if integer 
values are assigned to £ and n and increments A£ and An are chosen to 
unity. This gives rise to uniform spacing in the transformed plane, 
with no loss of generality since these increments cancel from the 
equation anyway. 

The generating system of equations (2.5 - 2.6) is represented by 
second-order central finite difference approximations in the transformed 
plane. Thu quantities A£ and An disappear by cancellation in all 
difference equations. Equations (2.5) and (2.6) are solved by the point 
successive over-relaxation (SOR) , [33], scheme after control functions 
P and Q have been specified. 

2 . 3 Control Functions 

The inhomogenous terms (P & Q) in equations (2.3) and (2.4) can 
be automatically chosen to obtain control of spacing, orthogonality 
and stretching. 

2.3.1 Th ompson et.al. Approach 

Thompson's approach consists of determining a correspondance between 
n values and the radius of concentric circles distributed between two 
circles with radius r^ and r 2 , one circumscribing the airfoil and the 
other tangential to the outer boundary respectively. Applying the coor- 
dinate generating equations (2.5 - 2.6) to the radii r^ and r 2 , while 
noting that n “ 1 on the airfoil and n a JL on the outer boundary, 
results in the following expression for Q: 

12 





Q(S,n) 


( 2 . 11 ) 


r r"(n) _ r ' (o) . 
l r'(n) " r (n) 

where r(n) is a function of the hyperbolic tangent [34], The effect of 
Q is to place a line corresponding to n *• k at a distance proportional 
to the r^ - rj from the body surface. Usually, to ensure proper 
resolution of the boundary layer the first line away from the boundary 
is placed at an approximate distance of one percent of the Blassius flat 
plate boundary layer thickness from the body, i.e. 

r(n=>2) - r, = 0.01( — - — ) (2.12) 

1 

In general, for the above approach, the control functions are determined 
from specified line distribution and the control functions have direct 
control over line spacing in the field (fig. 2a) . 

2.3.2 Sorens o n's Ap p roach ; 

Sorenson [35] determines the inhomogenous terms P and Q to control 
the spacing between mesh points, along mesh lines intersection the 
boundaries and the angles with which mesh lines intersect the boundaries 
(fig. 2b). P and Q are defined in terms of four new variables. In 
particular, for 1 < n < 1 max they are 

P(e,n> - p(5)e" a(n-1) + r(e)e ^ (2.13) 

Q(?,n) - q(5)e _b(n " 1) + s(£)e ^ (2.14) 

where a,b,c and d are positive constants and in particular a,b,c and d 
were set equal to 0.55. Note that terms p,q,r and s appearing in the 
above equations are functions of £ only. The four more unknowns p,q,r 
and s Introduced in equations (2.5) and (2.6) through equations (2.13) 


13 







and (2.14) require four equations. At the inner boundary ( ra 1) , the 
coefficient of r and s in equations (2.13) and (2.14) becomes very small 
and hence the wnd terms on the RHS of these equations can be dropped so 
that 


P(5,l) - p(0 

(2.15) 

Q(S,1) - q(0 

(2.16) 

boundary (n « n ) 

J max 

P(? * n max } “ r( ° 

(2.17) 

QU * r w ) " 3(0 

(2.18) 


Substituting equations (2.15) - 2.18) in equations (2.5 - 2.6) we 
obtain 


R,y - R.x 
r 1 n 2 m 
p(0 . { - 1 ] 


n=i 


,<o ■ 


ry*l 


(2.19) 


( 2 . 20 ) 


R.y - R,x 

t(E) . ^SjL—iJL, 


n°n 


max 


( 2 . 21 ) 


it 

‘I 

a 

*1 

R 


-R 3 y + R,x 

8(0 » 


n»n, 


max 


( 2 . 22 ) 


!• 

£ 
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where 


[ - (ax « 

- 26 % 

+ YK n-, )l n.l 

(2.23) 


- 2Sy in 

+ T! 'no ,1 n.l 

(2.24) 

(-(cix 

- 28 x 

+ YX ) ] 

(2.25) 

U 

Sn 

nn n=n 

max 

[ - (ay « 

- 2Py 

+ Y y„_) L „ 
nn n«n 

max 

( 2 . 26 ) 


Equations (2.19 - 2. 22) involve the derivatives at the inner and outer 
boundaries. At this point, if we assume that information about all 
these derivatives at the boundaries is readily available, we can compute 
the control functions P and Q using equations (2.13 - 2.14 ) for given 
values of £ and n in the field. 

The geometric constraints imposed by Sorenson will be used to 
define values of some derivatives at the inner and outer boundaries. 

The first requirement is that the spacing along 5 = constant lines 
between an inner boundary node at n = 1 and the corresponding next 
grid node at n “ 2 is specified by the user. Let this desired spacing 
in the physical plane be denoted by Asj^^ so' that we have 


As l n »i “ [<Ax) + (Ay) ] n»l 

in the limit Ax and Ay approach zero 


(2.27) 


ds 


ryl 


[(dx) 2 + (dy) 2 ]^ 


rpl 


(2.28) 


and transforming using the chain rule 

d3 ln“l “ + x n dn ^ 2 + + y n dn ^n»l (2-29) 


15 







for small distance ds along £ “ constant 


ds I = [(x 2 + y 2 ) z dn] , 

1 n^l n n it*! 


(2.30) 


or 


I r 2 2^ 

si , * (x + y J . 
»l 1 n=l n n n=»l 


(2.31) 


The second requirement is that the angle 0 of the intersection 
between the inner boundary and the £ = constant line is specified by 
the user. By using the definition of the dot product 


[ac • v n ] . = [|a?| |v n | cose] 


(2.32) 


or 


(5 n + £ n ] . = [ (£ 2 + £ Z )*(n^ + n^) “cose] t 

1 x x y y n°l x y x y J n«l 


(2.33) 


Using relations given in Appendix A. and equation (2.31) in (2.33) 
after some algebra we obtain 


n=l 


s (-x.cosO - y_sin0) 
r_Q_ — I _J l 

1 , 2 , 2 'it J . 

(x + y^) H"1 


(2.34) 


1 n»i 


s (-y cos9 + x sin0) 

r_IL — h S 1 

/ 2 . 2.H , 

(x + y ) n*i 


(2.35) 


Thus, the values of derivatives x^ and y^ at the inner boundary 
can be fixed by the user. Similar expressions can be obtained to 
fix the values of derivatives x and y at the outer boundary. 

After the values of x and y at the boundaries are specified, an 
iterative solution of grid generating equations (2.5 - 2.6) require 
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computation of the forcing functions P and Q in the field, which in 

turn, require information about derivative x_, x , x_„, x_ , x , 

* s V n’ ££’ 5n* nn* 

y<r» y . y»-rj y_ and y at the inner and outer boundaries. The 

'£ -V '55 'Sn 'nn 

desired values of spacing and angle at the boundaries, supplied as 
an input, will fix the values of x and y as discussed before. Also 
at the inner and outer boundaries q 15 constant, and hence values of 
x„, y„, x,.„ and y_,. are fixed. Derivatives x„ and y„ can be com- 

5 " V SS £n J Zn 

puted by differencing x and y with respect to £ and are fixed at all 

iteration levels. However computation of x and y at the inner and 

nn nn 

outer boundaries will require the use of one sided differencing schemes. 
These one sided finite difference approximations will require informa- 
tion about x and y at more than one point off the boundaries. As the 
values of x and y in the field will change with every iteration, the 
only derivatives that change with it are x and y . Thus, at each 
iteration level, the control functions P and Q change through x 
and y . The control over mesh spacing and angles in the field, 
introduced by equations (2.13) and (2.14), decays with the increase 
in values of (n-1) and (n -n) • The four control functions a,b,c and 
d in these equations determine the rate of exponential decay. It is 
interesting to note here that Sorenson's method is not overspecified. 
Since, control functions are to be determined we can specify additional 
boundary conditions. An iterative method such as SOR can be used to 
solve the system of governing equations. 

In the present investigation grids generated using Thompson's 
and Sorenson's techniques were employed. 
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Chapter III 

THE EQUATIONS OF MOTION 

The governing equations for a high Reynolds number incompressible 
flow field are the conservation equations for momentum and mass known 
as the Navier-Stokes equations. The governing equations in the present 
study are the time-dependent, incompressible, two-dimensional, Reynolds- 
averaged Navier-Stokes equations formulated in terms of the primitive 
variables. The pressure equation solved is the Poisson equation, 
derived by taking the divergence of the momentum equations. The Eulerian 
method is usually employed in computational fluid dynamics. This method 
involves a fixed control volume that is specified relative to a given 
coordinate system. Properties of the fluid are then specified as 
functions of both space and time. The conservation equations are 
approached using this methodology. 

3 . 1 Cons e rvation of Momentum 

For a given system, Newton's Second Law states that the rate of 
change of momentum is equal to the sum of the external forces acting 
on it. For an arbitrary material volume V, this law can be written as: 


III 


37 (pu i )dv + 


(Pu i )Ujn_jds 


<E « ‘ 1>6 ij ) “j da + a 


Pg ± dV 


( 3 . 1 ) 


The index "i" denotes any of the three cartesian coordinate directions 
x^, x^, x^, and the Einstein summation convention has been used for 
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I,,,*, -<r 


the index "j". The dimensional variables are: 


P “ density 
a velocity 

s « material surface 


n^ = unit vector, normal to s 
Z. , « shear stress tensor 

ij 

P = pressure 

= Kroncker delta 


** body-force acceleration 


The divergence theorem transforms Eq. (3.1) to 


ill 


5F (Pu i )dV + 


3^ (pu i u j )dV 


111 


E « - T?l + «i >d5 


Since this equation is valid for any arbitrary volume V 
integrands are continuous, the equation is 


A (P'u.) + TT“ (P«,u,) = 


JE- 


+ Pg, 


at a Xj i y ax.. ij ax t " & i 


For an incompressible flow density P is constant, so Eq. (3. 

3u i . a 


at 


, 3 / \ 1/ 3 s* 3p N 

+ (u i u j ) " p^;" z ij - a^ + g i 


(3.2) 

; when the 

(3.3) 

3) becomes 

(3.4) 
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and Stokes' hypothesis gives the shear stress as 

Z, 


ij 


3u 8u 

^ + 3 ^> 


(3.5) 


where P is the viscosity of the fluid. 

Equations (3.4) and (3.5) are normalized by defining the 
following dimensionless quantities. 


u . ® 
1 

u./U 
i » 

(3.6) 

V 

x^/Z 

(3.7) 

t = 

tu / z 
00 

(3.8) 

A 

P = 


(3.9) 

A 

P a 

p/o>i6 

(3.10) 

A 

8 l“ 

g ± Z 

■ 

(3.11) 


The reference quantities are 

U. ■ freestream velocity 

Z » characteristic length (airfoil chord 
for this case) 

P^ ** freestream viscosity 

Substituting equations (3.6) - (3.11) into equation (3.4) yields the 
normalized, time dependent, incompressible Navier-Stokes equations 

~ r + -r- (u ± u.) + ~ [p(— ~ + ~r^) ] + g t 

at ax. J ax. ax, ax, ax, 

j i J J i 

(3.12) 
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where Re is Reynolds number given by 


i*rutrt3rY'fi«fr «!■* »<•* 


Re - pUJt/p^ 


(3.13) 


Usually, viscosity is taken as constant for incompressible, non- 
conducting flow. However, in this study it is retained as a variable 
to facilitate the implementation of an algebraic model for turbulence. 

3 . 2 Con serv ation of Mass 

For a. given system in which matter is neither created or destroyed 
the law of mass conservation (continuity) can be written as 


' f 


i£ 

at 


dv + 


rr 

JJ 


pu.n ds 

J J 


(3.14) 


Applying the divergence theorem and eliminating the volume integrals 
as before, Eq. (3.14) reduces to 




which for incompressible fluids is 

3u . 

D = 


(3.15) 




(3.16) 


where D i3 the divergence of the velocity vector. Equation (3.16) 
remains unchanged by the introduction of non-dimensional variables. 

3.3 The Pressure Equation 

The irtcompressiblity constraint eliminates the equation of state, 
which relates pressure, density and temperature. Hence, the real 
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difficulty in the calculation of the velocity field for incompressible 
flow lies in the unknown pressure field. The pressure, gradient forms a 
part of the source term of the momentum Eq. (3.12). Yet there is no 
obvious equation for obtaining pressure. The pressure field is indirectly 
specified via the continuity equation. When the correct pressure field 
is substituted into the momentum equations, the resulting velocity field 
satisfies the continuity equation. 

To obtain a Poisson equation for pressure a divergence operation 
is performed on the momentum Eq. (3.12). 


a 2 P_ /“i /Vd , 3 \ 

2 dt3x 3x.3x. 3x.3x. 

3x i i j i j 


+ P 



3t3x ± 


(3.17) 


Substituting Eqs. (3.5) and (3.16) into (3.17) leads to 


„ 3u. 3u 2 3u 3u 

2 .3D . i i . 3 u , i . i N 

p 3t M 3x, 3x. 3x.3x. '‘Bx. 3x/ 

i 3 1 J 1 i 


+ 2 


^7’S 


(3.18) 


3 D 

where it is assumed that D = 0, but — $ 0 and that 

0 t 


3x . 


= 0, i.e. , 


the body-force acceleration is applied uniformly to the entire field. 

In deriving Eq. (3.18) D can be extracted and set equal to zero 
and thus will be zero; however, due to computer round-off error ~ ■ 

d t d C 

is expected to retain an appreciable value. Therefore the derivative 
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— serves as a corrective term to adjust the pressure in an 
at 

effort to satisfy the. continuity equation, as suggested by Hirt and 
Harlow [36]. 

з. 4 Normalized Governing Equations in Two Dimensions 

From this point on, all variables used will be non-dimensional, 
and the circumflex (~J will be dropped from the notation. Carrying 
out the indicated summations and identifying u. , u^, x^ and X 2 with 

и, v, x and y respectively yields the two-dimensional governing equa- 
tions. 

It should be noted the momentum Eq. (3.12) is written in conserva- 
tive form. As shown by Roache[9], this conservative form allows the 
finite-difference equations to preserve the Gauss divergence property 
of the continuum equations. Also, the Rankine-Hugonoit shock relations 
were derived using the conservative form. Thus, shock jump conditions 
are automatically satisfied since the conservative variables are 
continuous across the shock and need no special treatment because of 
discoritinuitites. Since the flow under investigation in this research 
contains no such discontinuities a., further simplification can be 
obtained using the non-conservative form. The non-conservative form 
of Eq. (3.12) is obtained by expanding the convective derivatives 
and using the continuity equation (3.16). 

(u 2 ) + (uv) = uu + vu (3.19) 

x y x y 

(v 2 )^ + (uv) x = Wy + uv x (3.20) 

Thus, in cartesian notation the governing equations in non-conservative 
form and two-dimensions are 
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u t + uu x + vu y “ - p x + H tlJv2u + 2 Vx + M y (u y + v x } 3 + h (3>21) 


1 


t + UV X + W y = - p y + Ri- t,J7 v + 2U y V y + M X (u y + V 3 + 8 2 (3 ’ 22) 


y y 


V 2 p = -D - (u 2 + 2u v + v 2 ) + — [y V 2 u + y V 2 v 
r t X y x y Re x y 


+ y u + y. (u + v ) + y v] 
xx x *xy y x yy y 


(3.23) 


8 • ^ The T ran sformed Equation s 

Equations (3.21), (3.22) and (3.23) in the physical xy plane are 
transformed into the 5n-plane using the definitions and relations 
given in Appendix A. The individual components of the transformed 
equations which are valid on a rectangular field (or a combination 
of rectangular fields) in the Cn-plane are 


[u t ] = [u t J - (u x x t + u y t ) = [u t ] (3.24) 

x,y ?,n y S,n 


[v t ) = [v t ] - (v x x t + v v t ) = [v t ] 


x,y ?,n 


€,n 


(3.25) 


u x • ( v? - y? u n )/J 


(3.26) 


u * (x,_u - X U_)/J 

y 5 n n v 


(3.27) 


v x - ( Vc - h\ )/J 


(3.28) 


v y " (x £ v n ‘ x nV /J 


(3.29) 


p ■ (y p„ - vn ) /J 
*x / n l 4 ' 5 n 


(3.30) 
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P y = ( Vn “ X n p 5 )/J 


? 2 
V u - (ctu__ - 2gu_ + yu + ou + xu_)/J 
55 5n nn n 5 


2 2 
V v = (av_ c - 2$v + yv + ov + tv_)/J 

55 5n nn n 5 


u x = (y n u 5 " V 5 y n )/J 


= / J 

y 5 n n 5 


ts i ] = [g i ] P “ (g i x t + g i V = 

x,y 5,n x y 5,n 


[g 2 ^ = Is 2 l - (s 2 x t + §2 y t^ = ^ s 2^ 

X,y 5,n X y 5,n 


» 2 p - <«P,, - 28p ti) + YP nn + »P n + tp { )/J 2 


[B t ] ■ [Dl - (VS; + D y > - [DJ 
x,y 5,n y 5,n 

2 2 2 
U = [y \i~r- ~ 2y v u + y u - (y y __ - 2y^y y r 

xx y n 55 5 n 5n 5 nn n 55 C n 5n 

+ y^y )n - (y^x._ - 2y y x + y^x )u ]/J 3 

5 nn y y n 55 5^n 5n *5 nn x J ' 


P = {(x,.y + x y _)u_ - x_y^p - x y u._ 

xy 5 n n 5 5n 5^5 nn n y n 55 


(3.31) 

(3.32) 

(3.33) 

(3.34) 
' 3 . 35 ) 

(3.36) 

(3.37) 

(3.38) 

(3.39) 

(3.40) 


+ [x y x r _ + x y x - (x y_ + x.y )x,. ]p 

n y n 55 55 nn n y 5 5 n 5 n x 


+ [ Vn y 55 + X 5 y 5 y nn " (x n y 5 + X 5 y n )y 5n ]y y }/J (3 ‘ 41) 
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(3.42) 


2 ? 2 
H » [x U„,. - 2x r X M, + x"M - (x y_. - 2x_x y. 
yy n 55 5 n 5n 5 nn n*55 5 n'5n 

2 2 2 2 
+ x„y )p - (x x. r - 2x.x x r + x r x )p ]/J 
5 nn y n 55 5 n 5n 5 nn x 4 

where a, 3, y and J are defined in equations (2.7) - (2.10) and a and 
T are given by 

a « {-y (aPx^ + yQx^) + x^aPy^ + yQy^)]/J (3.43) 


t = [y^aPx '+ yQx^) - x^(aPy + yQy^l/J 


(3.44) 


The discretization of the transformed versions of equations (3.21), 
(3.22) and (3.23) and the numerical procedures used to obtain their 
solutions are discussed in Chapter V. 

3.6 Turbulence Model 

Since the flow fields of interest are turbulent, the solution of 
the Navier-Stokes equations must take into account the effects of the 
random fluctuations of the dependent variables inherent to turbulent 
flows. The turbulent nature characteristics of these flows can be 
accounted for in the numerical soltuion by a variety of eddy viscosity 
models ranging from locally dependent algebraic models to the more 
complex higher order closure models. A paper by Marvin [19] provides 
a comprehensive survey of turbulence models generally employed in 
computation of external aerodynamics flows of practical interests. To 
date no single turbulence model has emerged that can be applied to the 
variety of flows encountered in computational aerodynamics. Also, the 
use of higher closure models will not necessarily give more accurate 
solution. Therefore 't was decided to use the locally dependent eddy 
viscosity model. The turbulence model used in this research is an 
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extension of the Cebeci-algebraic viscoisty model [37] as modified 
and reported by Baldwin-Lomax [32 j. In this model distribution of 
vorticity is used to determine length scale which eliminates the 
somewhat uncertain process of finding the outer edge of the shear 
layer. The non-dimensional molecular coefficient of viscosity P in 
the laminar Navier-Stokes equation i3 replaced by 

|j *> 1 + e (3,45) 


f 

i 

! 

$ 


where e is eddy viscosity. The boundary layer region on a body 
consists of two layers, the inner layer and outer layer. The inner 
layer of this model accounts primarily for the laminar sublayer 
adajcent to the wall, with the outer layer accounting for the remainder 
of the boundary layer region. In cartesian coordinates the expression 
for the modified inner model based on Prandtl’s mixing length theory 
can be written as 


^ inner 



(3.46) 


where u» is defined as vorticity 


M 


3u 3v 
3y 3x 


(3.47) 


The mixing length in this model is obtained from Van Driest's sublayer 
model, and is given as _ - — 

y ww 

l s 0. 4y [1 - e 26 ^ ] (3.48) 

where y is the normal distance from the wall. 

The outer region eddy viscositymodel consists of a modified 
closure-type model defined by the equation 


i 
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(£) outer “ °- 0268p F l F 2 ( y } 


(3.49) 


where ? 2 (y) is the Klebanoff intermittency factor given by 

F-(y) - [?. + 3.5 (-it^) 6 ] 1 (3.50) 


max 


and 


F, = y F 
1 max max 


-y /FT 


w w . 


F(y) = y ! iu | [1 - e 


26 p 


(3.51) 


(3.52) 


The quantity F is the maximum value of F(y) that occurs in a 
’ max 

profile and y is the value of y at which it occurs. 
r 'max 

The eddy viscosity in the wake region is given by the equation 
(3.49) with F^ and F(y) defined as 


F, = 0.25y u /F 
1 'max dif max 


F(y) * y|wj 


(3.53) 

(3.54) 


where 


, _ 2" _ 2. 

u dif " ( u + V ^max - ( u + v ) min 


(3.55) 


For some cases under investigation, the boundary layer transition 
points were set by assuming that transition occurs at the minimum pressure 
points and for the other cases, the transition points were maximum air- 
foil thickness points on the airfoil surface. The above prescribed twor- 
layer eddy-viscosity model was successfully used by Baldwin and Lomax [32] 
and other investigators to predict separated flows. 
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Chapter IV 


BOUNDARY AND INITIAL CONDITIONS 

Boundary and initial conditions must be defined in order to solve 
the governing partial differential equations of a given flow field. 

Since important features such as boundary layer arises from boundary 
conditions, these conditions must be carefully defined. The conserva- 
tion equations for incompressible flow about an airfoil when formulated 
in terms of primitive variables require initial velocity and pressure 
distributions and either Neumann or Dirichlet boundary conditions for 
the velocity and pressure on the boundaries. 

4 . 1 Init i al Conditions 

Since the governing equations contain time dependent terms, 
initial conditions must be specified for the solution to proceed. 

Initial values of velocities and pressure must be imposed over the 
field. The values of non-dimensional velocities and pressure were 
set to zero at a time t «* 0. Once an initial case for the flowfield 
had numerically converged to a valid solution, each succeeding time 
step was initialized by using velocity and pressure distribution of 
the preceding time steps. 

4 . 2 Free-stream and Downstream Boun d ary Conditions 
Computationally, the free-stream boundaries are generally placed 

at a reasonable distance from a body such that uniform flow conditions 
remain undisturbed by the presence of the body. Velocities and pressure 
are completely specified at the free-stream boundary (Section ^ , 
in Fig. 1),. 

The flow Is accelerated from zero to a desired final velocity 
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using the body force terms and g^ in equations (3.21) and (3.22). 
The values of the velocities on the free-stream boundary during the 
acceleration phase were determined in the following manner 


gjdt 

0 < t < 1 

(4.1) 

s 2 dt 

0 < t < 1 

(4.2) 



(4.3) 


where = cos 1 p and g ^ - sin \p, where ip is the angle of attack. 

For each time step of the acceleration phase, the velocities on the 
free-stream boundary are found using eqs. (4.1) - (4.2) and are held 
fixed for computation of solution for that time step. After the 
acceleration phase the free-stream velocities are 

u^ = cos <p (4.4) 

v = sin tp 

The body force terms g^ and g 2 were set equal to zero after the 
acceleration phase. 

The boundary (T^) is placed a great distance downstream. For 
this case, no velocity gradient exist at the downstream boundary. 

The pressure at the downstream boundary was set equal to the free- 
stream pressure. These boundary conditions can be written as 

u • 0 (4.6) 

x 

v n 0 (4.7) 

x 

p - 0 (4.8) 
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Also the effects of different forms of downstream boundary conditions 
investigated will be presented in Chapter VII. 

4 . 3 Body-Surface Boundary-Conditions 

The airfoil surface is considered to be a no-slip and impermeable 
boundary. The no-slip and no-transp:Lration conditions at the airfoil 
surface can be written as 

u = 0 (4.9) 

and v = 0 (4.10) 

The pressure on the airfoil surface is unknown, but can be approxi- 
mated using the normal pressure derivative in the following way. The 
momentum equations (3.21) and (3.22) are utilized to evaluate the 
normal derivative of the pressure. Due to no-slip no transpiration 
boundary condition at the surface, the transient and convective terms 
in the momentum equations drop out and we obtain 

= n • Vp = n • (g + —• V^u) (4.11) 

3n " Re 

where n is a unit normal and g is the body force vector. 

Initial attempts to use the above pressure boundary condition 

led to computational divergence. The simplified version obtained by 

* 

neglecting the viscous terms was used in the present study. 

n • Vp = n • g (4.12) 

The presence of the body force vector influences the pressure boundary 
condition during the acceleration phase; however after the acceleration 
phase is over equation (4.12) reduces to the familiar form 

n • Vp = 0 (4.13) 
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For the present case, the airfoil surface is represented by 
n ” constant line. The direction normal n can be given by Vq. 

Thus, equation (4.12) becomes 

Vq • Vp = Vq • g (4.14) 


or 


(£ n +5 n )P, + 

’X X y f. 


+ "X 


n x 8 l + n y g 2 


(4.15) 


Using Appendix A we obtain 

p n ” + J(_y f; 8 l + x ? s 2^ (4 .i6> 

Tlie surface pressure can be evaluated using a one-sided finite- 
difference approximation for P^. For the problem under consideration 
the boundary conditions at the airfoil surfaces are probably the most 
crucial. 

4 . 4 Re-ent r ant Boundary Conditions 

The re-entrant sections, and in figure 1 are not boundaries 
in the physical plane but represent points within the flow field. The 
branch cut is made between the trailing edge of the airfoil and the 
downstream boundary to eliminate discontinuity in the inner boundary in 
the transformed plane. The values of flow variables cannot be fixed 
at these boundaries but they should evolve as a part of the field 
solution. This insures the continuity in flow variables and their 
gradients across the cut. 

4 . .5 Trailing-edge Boundary Condition s 

In the transformed plane body surface is a continuous line; 
however^ in the physical plane the trailing edge is a sharp point. The 
surface-normal vector Vq is discontinuous at the trailing edge. This 
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geometric discontinuity leads to unequal trailing edge pressure 
found using equation (4.16). The basic assumption that there be 
no unbalanced forces at the trailing edge would be violated. To avoid 
this problem, the trailing edge pressure was found by taking the average 
of the trailing edge pressure on the upper and lower airfoil surfaces 
(points (NWE,1) and (NWS,1) in Fig. 1). However it led to jump in the 
pressure at. the trailing edge which is physically unrealistic phenomena. 
To obtain smooth pressure distribution at the trailing edge, the follow- 
ing extrapolates were found useful in the present study 

P(NSW,1) = ^(P (NWS +1,1)+ P(NWE - 1, 1)) (4.17) 

P(NWE,1) = %(P(NWS +1,1)+ P(NWE - 1, 1)) (4.18) 
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Chapter V 




SOLUTION ALGORITHM 

5 . 1 Numerical Procedure 

The governing equations are the two dimensional, time dependent 
Navier-Stok.es equations in the non-conservative form. The Poisson 
equation for the pressure is obtained by taking the divergence of the 
momentum equations and utilizing the continuity equation. The two 
momentum equations and the Poisson pressure equation form a set of 
three, governing equations for three flow field unknowns u, vand p. 

These governing equations are solved in the transformed plane 
for each field node using a fully implicit finite-differencing algo- 
rithm. This implicit algorithm is obtained by means of backward- 
time and central-space differencing of derivatives in the transformed 
plane. The governing finite-difference equations in an implicit 
form are fully vectorized and solved simultaneously at each time step 
using a checkerboard matrix iterative technique (Chapter VI) . 

5 . 2 Finite-difference Approximations to Governing Equations 

As discussed before in Section 3.5 the task of obtaining the 
transformed governing equations is straight forward and requires 
substitutions of the transformed expression for the derivatives (3.24 - 
3.44) in the governing equations (3.21 - 3.23). The presentation of 
fully transformed governing equations has been avoided here for simpli- 
city; however, this section will detail the specifc transformations that 
are pertinent to the final form of computational equations. 

All spatial derivatives in the transformed equations are approxi- 
mated by second-order-accurate central-difference ixpressions as follows: 
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M „ fi+LJJllizLi 

3C 1§J 2A 5 


3 2 f „ f i+l,,1 " 2£ i,J + f i-l,j 

2 2 
35 i.j A5 


( 5 . 1 ) 


(5.2) 


Similar finite-difference expressions are used to approximate n 
derivaitve. The second-order accurate expression for the cross- 
derivative is 


JUL , fj j U±i ~ ii+LtizLl * £ i-l 

3£3n 4A5&0 


(5.3) 


The grid spacing A£ and Ari is chosen to be unity because of the construc- 
tion of the mapping from the physical plane to the transformed plane. 

As presented in section 4.2, the flow is accelerated from rest 
to the final desired free stream velocity. Hence, the temporal 
derivatives are represented by the first-order-accurate two-point 
backward-difference scheme at the first time step and by second-order 
accurate three-point backward-difference at all subsequent time steps. 

The expression for two-point backward difference is 

(n) 

i»J 


3f 

3t 


.(n) (n-1) 

At 


(5.4) 


and for the three-point backward difference is 


ii 

3t 


(n) 3f5 n) . - 4f (n “ 1) + fj n ' 2) 


(5.5) 


where the superscript (n) indicates the time level. 
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To obtain the computational form of the governing difference 
equations for an iterative scheme, we must combine the diagonal terms 
(those with subscript i,j) of spatial derivatives with an appropriate 
temporal derivative term. As the central-difference approximations have 

been used for the spatial derivatives, the terms with subscript i,j 

2 

will appear only due to tranformation of V ( ) in equations (3.21 - 

2 

3.23). For completeness, we transform V u in the u momentum. From 
equation (3.32) or Appendix A we have 
2 1 

V U = — r- (oui . - 2Bu_ + YU + OU + TU„) 

j2 ££ 5n nn n £ 

The finite-difference approximation of derivatives u„ and u will 

££ nn 

involve the diagonal terms and approximation of all other derivatives 
ill involve off-diagonal terms. Separation of diagonal and off 
diagonal terms gives 

V 2 u = C V 2 u) D + (V 2 u) qd (5.6) 

where 

(V 2 u) d = - — (a + y)u i j (5.7) 

J 

(V 2 u) od = ^ ta(u i+1J + u^) + Y(u iJ+1 + u^j^) 

-2$u r + au + Tu r ] (5.8) 

£n n £ 

Note that in the remainder of this section terms such as u , v , etc. will 

x’ y ’ 

appear but are to be implicitly assumed to have been evaluated according 
to equations (3.24 - 3.44) or relations in Appendix A. Substituting 
equation (5.6) in the u momentum equation (3.21) and combining the 
diagonal term with the temperal term at time level n, we obtain the 
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computational form of the difference equation. 


[ 4 - + - (a + Y)]uj ll) = B + (RHS) 

/Re X » J ' X 


(5.9) 


where 


(RHS), ° -uu - vu - p + — [M(V 2 u)_ n 
1 x y *x Re OD 


and 


+ 2P u +h(u + v)] + g, 
XX y y x' i 5 1 


A = 1 


B 


(n-1) 

At 


for two-point backward differencing, or 

3 


B 


(n-1) (n~2) 

4u - u, . 

*.J 

2At 


(5.10) 

(5.11) 

(5.12) 


(5.14) 


for three-point backward differencing. The computational difference 
equation for the v momentum, derived in similar fashion is 


[- 4 - + -jr~ (« + Y)]v (n) = C + (RHS), 

/Re 


(5.15) 


where 


(RHS) 2 = -uv x ~ w y - p y + ~ [m(V 2 v) O0 
+ 2 u y v y + U X (u y + v x )] + «2 


(5.16) 


and 


or 


(n-1) 

V i 1 

C = — — for two-point backward differencing 


4v (n_1) - v (n_2) 

Z1A J.J— 

2At 


(5.17) 


for three-point backward differencing 

(5.18) 
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The computational form of the Poisson pressure equation can be 


obtained in a similar manner from equation (3.23). The term D t 

represents the time derivative of the divergence of the velocity 

vector. It is assumed that the conservation of mass is satisfied 

at: the most: recent time level (i.e., D n = 0); however values of the 

divergence at previous time levels have been retained as a corrective 

term. Thus 

D (n-D 

D t = - — for two-point backward differencing (5.19) 

-4D^ n_1 ^ + D^ n ~ 2 ^ 

D t = 2 ^ for three-point backward differencing 

(5.20) 

and the computational pressure equation takes the following form 


J lp ' ■ ( ’ 2 P)oo + \ + “* + V, 


2 2 2 2 
+ v - — p V u + M v 
y Re x y 


+ Pu+U (u + v ) + M v ] 
xx x xy y x yy y 


(5.21) 


where D t can be approximated using either equation (5.19) or equation 
(5.20). The first approximation is first order accurate while the 
second approximation is second order accurate. Several computer runs 
were made with the two approximations for comparison. No significant 
difference were found between results obtained using the first order 
and second order accurate approximations. Also, for particular test 
runs, none of these two approximations was specifically responsible 
for decay or divergence of solution. The above tests were not entirely 
conclusive. 
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The goal is to find the steady state solution regardless of 
accuracy of the transient solution. Since the time derivative terms 
will hopefully disappear in the steady state and higher order approxi- 
mations usually require more operations per mesh point, we used first- 
order two-point approximation for D t in the pressure equation (5.21). 

5 . 3 Fin it e-difference Approximations to Boundary Conditions 

The downstream boundary condition equations (4.6 - 4.7) are 
transformed according to the relations in Appendix A. Equation (4.6) 
for the lower downstream boundary (£=1) takes the following form 


u = — {y u_ - y,.u } = 0 


u f = (-% u 

e y n n 


(5.22) 

(5.23) 


Using one sided three point forward differencing for u 


1 y E 

u i , = [4u 0 . - u_ , - 2(— )u ] 

l,j 3 2, j 3,j > ' n 


(5.24) 


Similarly fc the upper-downstream boundary condition (£ = IL) using 
three point backward differencing for we obtain 


‘iLJ - 5 ' u IL-2,j + 


(5.25) 


[ 

* 


i 

i 



In the above equations, derivative is evaluated using central- 
difference approximations. Replacing velocity u by velocity v in 
equations (5.24) and (5.25) we can obtain expressions for equation 
(4.7). 

Pressure values on the airfoil surface are determined using the 
Neumann boundary condition (4.16). Using a three point forward- 
difference approximation for p the expression for airfoil pressure is 




v & -A ,4-. -*<K 1-^ $£ .f ^ 
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(5.26) 


5 . 4 The Re-entrant Boundary 

The procedure for the evaluation of flow field variables (u, v 
and p) on the cut extending from the outflow section (fig. 1) deserves 
special attention. The two re-entrant sections and r,. resulting 
from the cut are one and the same lino in the physical plane. Thus 
corresponding points on the two re-entrant sections have the same x,y 
coordinates in the physical plane but different £ values. The momentum 
and pressure equations on the re-entrant section can be solved assuming 
continuous derivatives across the cut. However, in this study, flow 
variables on the re-entrant section were found by averaging the 
corresponding values above and below the branch cut. As the grid 
spacing required at the branch cut to resolve the flow is very small 
in a C-type grid, the averaging gives almost the same values of the flow 
variables as those found solving the governing equations at the re- 
entrant section. An expression obtained using notations of fig. 1 for 
two nodes in the computational plane that correspond to the first node 
off the trailing edge in the physical plane is 


f NWS-l,l “ f NWE+1 , 1 “ ^ f NWZ+l,2 + f NWS-l,2^ (5.27) 

Similar relations were used to find the values of the velocities and 
pressure at nodes on the branch cut. This approach simplifies compu- 
tation of flow variables at the re-entrant section without sacrificing 
accuracy and hopefully enhances the computational efficiency due to 
less involved operations. 
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5 • 5 Ar tifici a l Viscosity 

Central-differencing schemes frequently display oscillations on a 
coarse grid. The present implicit scheme exhibits oscillatory behavior 
at high Reynolds number due to inaccuracies introduced by finite- 
differencing. Unless these extreme oscillations are damped out the 
numerical solution becomes useless. In most cases under investigation, 
the solution started diverging about time t = 1.0 without inclusion of 
artificial viscosity. The use of artificial diffusion was found necessary 
to obtain steady state solution. The pressure oscillations were 
responsible for fluctuations and discontinuities in the velocity field. 

One possible source of the pressure oscillations was the divergence of 
the velocity vector which is a part of the source term of Poisson pressure 
equation (3.2.3) and may have retained significant magnitude. The basic 
assumption to obtain the Poisson pressure equation was the preservation 
of the continuity at the most recent time level. Thus, significant 
deviation from the satisfaction of the continuity equation can contami- 
nate the pressure field. Incorporation of artificial viscosity based on 
the divergence of velocity at the current time level can damp extreme 
oscillations. The modified non-dimensional viscosity coefficient in 
the momentum equations (3.21 - 3.22) is given by 

U - 1 + e + M (5.28) 

a 

where e is eddy viscosity term discussed in section 2.6 and is 
artificial viscosity. One possible form of the artificial viscosity is 

U a “ QReJ | V • V| (5.29) 

Note that this form of artificial viscosity has units of eddy 
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viscosity and it has the advantage of being analytically zero. The 
term fi can limit the effects of artificial viscosity and will be a 
constant or a variable derived from the flow characteristic. This 
particular form of artificial viscosity has a desired property of 
being proportional to |v • v| and only becoming effective in regions 
where the divergence of velocity is significant. The fluid dynamics 
phenomena investigated with various form of artificial viscosity will 
be discussed in. Chapter VII. 

Strictly speaking, the computation of viscous derivatives for the 
momentum equations (3.21 - 3.22) and the pressure equation (3.23) should 
use the modified viscosity coefficient (eq. 5.28) . However this approach 
lead to divergence of the solution and hence the viscous derivatives 
for the governing equations were computed using the viscosity coefficent 
given by equation (3.45). The artificial viscosity u was incorporated 

a 

in the viscosity coefficient p of the momentum equations at every time 
step . 

5 . 6 Smoothers 

At early time stages, the solution may contain enough noise to 
excite oscillations. Nonlinear interaction will amplify these 
oscillations which in turn may destroy the solution. In such cases, 
we wish to filter out the unwanted oscillations from the solution. In 
most test runs, wavy divergence of velocity field was obtained in the 
direction of ( ■ constant lines. If somehow, a smooth divergence of 
velocity field can be obtained, It can reduce the pressure oscillations. 
Two types of smoothers, one using the adjacent nodes in 5 » constant 
direction and the other using the four neighboring nodes were investiga- 
ted. The latter gave better overall smoothing of divergence of velocity 
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field. The expression for the smoother is 


i,;i 


iff 

2 U i,j 


+ 4 (f i+l,j 


+ f 




+ f iJ+l + 


ij-l 


)] (5.30) 


On the other hand, pressure exhibited excessive oscillatory be- 
havior in the direction of n ■ constant lines. Smoothing of the pressure 
itself, using equation (5.30), lead to Incorrect pressure solution, but 
smoothing of the source term of the Poisson pressure equation (3.23) can 
smooth out the pressure field oscillations. Denoting the source term by 
S, with the assumption that all terms on the right hand side of equation 
(3.23) are lumped into S, we can smooth out the source term by using 
S instead of f in equation (5.30). The results obtained using diver- 
gence of velocity and source term smoothers were almost the same. As 
computer operations for the source term smoother are more involved and 
smoothing operation is usually required at every iteration it was not 
investigated further in the light of computational efficiency. 

The divergence of velocity smoother, applied at every time-step, 
was employed to reduce, the need of: artificial viscosity. 
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Chapter VI 

VECTOR PROCESSORS ANI) CHECKERBOARD SOR 

In the last decade, significant progress has been made in the area 
of algorithms that are used for solving the governing flow field equa- 
tions. Computer codes employed in many engineering applications still 
use large amounts of computer resources. The basic requirement is that 
the algorithm be efficient. In practical terms this means obtaining the 
solution with desired accuracy using the least amount of computer resources. 
Any improvement in the numerical scheme used can certainly enhance the 
efficiency. Some improved algorithms have been mentioned in Chapter 1 
with appropriate references. Furthermore, in many cases the computers 
available play an important role in the development of efficient algorithms. 
Hence the computer achitecture such as serial, vector or parallel certainly 
dictate the basic requirement of algorithms. Frequently, the structure 
and size of computer memory and data mangement system can play a crucial 
role in the implementation of efficient algorithms. 

6 . 1 Vector Processors 

The advent of high performance sixth generation computers such as 
the CYBER-200 and CRAY-1 series, provides an important breakthrough for 
computationally demanding engineering problems. These supercomputers 
incorporate vector processing capabilities to provide the computational 
power required by large scale numerical simulation [21]. 

Vector processors are generally divided into two main classes of 
architecture: memory to memory (MM) and register-to-register (RR) . 

Normally MM architecture vector processors operate at its highest level 
of performance when algorithms being processed have the following 
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characteristics. Operand and result vectors are stored contiguously 
in memory i.e., successive elements of the vector must be stored' in 
adjacent memory locations. The length of the vector is long. The 
example of: MM architecture is the CYBER-200 series machines. RR 
achitecture usually involves some type of cache between main memory 
and functional units. The fundamental idea of cache, organization is 
that by keeping most frequently accessed instructions and data in the 
fast cache memory, the cache is only a small fraction of the size of 
main memory. If the active portions of the program and data are placed 
in a fast cache memory, the average memory access time can be reduced 
considerably, thus reducing the total execution time of the program. 

The cache is the fastest component in the memory hierarchy and approaches 
Che speed of CPU components. These types of vector processors usually 
achieve their highest level of performance when processing algorithms 
that satisfy the following requirements. Parallel execution of the 
functional units is maximized. The example of RR architecture is the 
CRAY-1 series machines. 

The above-described two types of vector processors are called 
pipeline processors. Pipeline is a technique of decomposing a sequen- 
tial process into subprocesses with each subprocess being executed in a 
special dedicated segment that operates concurrently with all other 
segments. A pipeline can be visualized as a collection of processing 
segments through which binary information flows. Each segment performs 
partial processing dictated by the way the task is partitioned. The 
result obtained from the computation in each segment is transferred to 
the next segment in the pipeline. The final result is obtained after 
the data have passed through all segments. The name "pipeline" implies 
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a flow of information analogous to an industrial assembly line. It 
is characteristic of pipelines that several computations can be in 
progress in distinct segments at the same time. 

The CRAY-1 series are RR type pipeline machines which operate most 
efficiently on vectors which are of length 64 or a multiple of 64. 

The reason is that the vector registers hold 64 words which are sent 
to the pipeline. Thus in the CRAY-1 series machines the vector registers 
are limited to 64 elements and hence extremely long vector lengths will 
not necessarily enhance the computational efficiency. The CRAY-1 
memory section normally consists of 16 banks of memory. The memory 
size can be as large as about 1 million words. Each word contains 44 
data bits and 8 check bits. The control of data flow between the parallel 
functional units and hierarchically organized memories is of significant 
importance for algorithm efficiency. 

The CYBER-200 series are MM type pipeline machines which operate 
more efficiently as the vector length increases. Each vector instruction 
involves a startup time, the time required to produce first result. 

Since startup time becomes relatively less important as the vector 
length increases, the vector operations become more efficient. Thus it 
is desirable to work with moderate to long vectors on the CYBER-200 
series machines. The CYBER-203 has about 1 million words of primary 
memory with virtual memory architecture. Memory on this machine is 
called as pages, which are of small and large size. The small page is 
made up of 521 words of 64 data bits and the large pages are of 65,536 
words. A user can have access to about seven large pages in primary 
memory at a time. The movement of data from secondary memory into 
primary manory . involves moving of pages. This movement of pages in 
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in and out of primary memory is called page fault and involves startup 
time and transmissio \ time. It is desirable to make most efficient use 
of data when it is in primary memory to avoid situations when the machine 
time spent on data management makes up a considerable part of the total 
time. 

Thus performance on a vector processor can vary widely as a function 
of algorithm, implementation and data management. 

6 . 2 Checkerboard SOR 

As the computers discussed in the above section attain their highest 
level of performance when processing vectors, it is clearly desirable to 
search for methods that can take advantage of the vector operation 
capabilities without suffering significant loss in convergence rate 
compared to widely accepted.methods for serial computers. 

Generally the choice of an appropriate algorithm is dictated by 
whether the flow is subsonic, transonic or supersonic. Although it is 
the steady state solution that is generally sought one often uses 
the time dependent equations to reach steady 3tate. An explicit 
algorithm which can be easily vectorized may have much slower convergence 
rate. With explicit methods entire two or tnree-dimensional grids can 
be considered as one long vector. On some machines this will lead to a 
high level of optimization. The solution of the three dimensional com- 
pressible Navier-Stokes equations obtained using vector processors have 
been published by several investigators. Smith et.al. [38] and Shang 
et. al. [39] solved these equations using an explicit scheme on the CDC 
STAR-100 and CRAY-1 computers respectively. For the 3-D problems solved 
using an explicit scheme, the vector lengths were restricted to the 
number of grid points in each 2-D plane due to efficient use of computer 
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architecture in Shang's investigation and due to efficient data manage- 
ment for memory in Smith's study. Spradley et. al. [40] solved these 
3-1) equations using general interpolants methods (GIM) on the CDC STAR- 
100. He chose weight functions suc-h that to produce explicit finite- 
difference type analog and used the vector lengths equal to the total 
number of grid points in a 3-D flow field. 

AJ.thou.gh long vectors available at each time step for explicit 
schemes may increase efficiency of some vector processors, the large 
number of time steps required to reach the steady state may adversely 
affect the overall performance of the, algorithm. Furthermore, in many 
cases, one is only interested in obtaining the steady state solution 
as fast as possible without regard to the accuracy of the transient 
solution. The time step restriction imposed by stability consideration 
is a major disadvantage of explicit' schemes. Hence there is increased 
interest in implicit schemes in recent years. Also, for implicit 
schemes, one frequently uses the time dependent equations, and fairly 
accurate steady state solution is reached with larger time steps. 

Although implicit methods are usually linearly unconditionally stable, 
however there exists time step restriction based on accuracy require- 
ments. Also operations for an iterative relaxation procedures are 
more involved. The development of an efficient relaxation method is an 
important element for implicit algorithms. 

The most widely used classical relaxation methods are generally not 
suitable for vector computers. The point accelerated successive relaxa- 
tion (SOR) method [33], which is perhaps most frequently used, is reli- 
able and very competitive for many problems. The convergence rate of point 
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SOR depends partly upon using updated values at adjacent points while 
solving for a given point. 

Point SOR schemes can be efficiently implemented on a scalar 
machine. However, for vector processors, vectors must be stored and must 
be available for concurrent computer functions required for desired 
arithmetic operations. This requirement is very restrictive and the 
classical point SOR method is not suitable for vector processing in its 
original form. There are some possible ways of system ordering for 
solution of PDE using a rectangular grid on a vector machine. Suffi- 
ciently large vectors can be identified within the field or subfield 
by (a) associating vectors with alternate rows or columns (ZEBRA) or 
(b) associating vectors with alternate field points (red - black) . 

Option (b) is a simple way of making point SOR suitable for vector 
processing. This modified SOR is usually referred to as checkerboard 
SOR or hopscotch method in case of parabolic problems. 

Early work related to the hopscotch method was presented by Cordon 
[41] in 1965, many years before vector processors became available. 

His work was motivated by the favorable stability properties of the 
method. Gordon(41] described the original technique as "A non-symmetric 
difference equation" obtained using explicit and implicit finite dif- j 
ference schemes at alternate mesh points and showed that combined scheme 
was unconditionally stable. Scala et. al. [42] applied this technique 
to solve the Navier-Stokes equations about a circular cylinder using 
a cylindrical. coordinate system on a serial computer. Gourlay, et. al., 
[43, 44] presented the original technique in a more general form and 
showed tnat the checkerboard method can be regarded as an Alternating 
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Direction Implicit (ADI) method with the coefficient matrix split 
in .a special way. The fundamental idea of an ADI scheme is of 
splitting the problem into a series of simpler problems. Normally, 
each simpler problem corresponds to each space dimension and in many space 
dimension problems complexity increases considerably. The major advan- 
tage of the checkerboard algorithm is that it can be always decomposed 
into two simpler problems (two stage process) irrespective of the number 
of space, dimensions. 

For illustrative purpose, it is convenient to consider a simple 
model problem. Some detail for solving the Poisson equation using 
checkerboard-SOR will be presented. Let us consider the Poisson equation 


2 2 

— J + — J = S(x,y) (6.1) 

3x 3y 

with simple Dirichlet boundary condition on the boundary and with the 
field subdivided in square cells to length h as shown in fig. 3. 

Using central finite-difference approximations at a mesh point (x^, y^) 
equation (6.1) can be written as 


4f 
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i+l>j 


- f 
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i-l>j ij+l ij-l 


-h 2 S 


i»j 


( 6 . 2 ) 


Instead of considering natural ordering of mesh point, i.e. sweeping 
rowwise, let us visualize the field mesh points as forming a red-black 
chess board. This red-black ordering can be defined as follows: Cell 

field mesh point (i,j) red if (i+j) is odd and point (i,j) black if 
(i+j) is even. Hence the red unknowns will be the set of all f^ ^ for 
which (i+j) is odd and similarly for the black unknowns when (i+j) is 
even. Applying the classical SOR to the red and black unknowns it can 
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be shown that each classical SOR iteration can be split into two stages. 
The first or red stage consists of improving the red unknowns according 
to 

(n+l,r) JCy (n,r) 4. <r^ n »h) . r(n,b) 1 r( n »6) 1 <:( n »b)) 
f i,j 4 ( “ h S i,j + f i+l,j + f i-l,j + f i,j+l + f i,j-l 

+ (1 - K)f (6.3) 

and during the following second or black stage the black unknowns are 
improved according to 

<n+l,b) k ( 2 (n,b) (n+l,r) (n+l,r) (n+l,r) (n+l,r) 

f i,j 4^~ h b i,j + £ i+l,J + £ i-l,J + r i,j+l r i,j-l ; 

+ (1 -K)f $ n ! b> (6.4) 

In equations (6.3) and (6.4) e is a relaxation parameter used 
to accelerate convergence, superscripts n,r and b denote iteration 
level, red and black nodes respectively. During the red stage all red 
iterates are updated with the help of the adjacent black iterates and 
conversely in this particular case. Each state Is inherently parallel 
in that all iterates of the same color can be updated simultaneously 
without changing those of the other color. Each term on the RHS of equa- 
tions (6.3) and (6.4) can be represented as a vector, assuming scalar 
terms such as h, is a vector of desired length with each element of the 
vector having the same value. Let us assume k to be constant for the 
present case. Thus equations (6.3) and (6.4) can take advantage of 
vector processing capabilities of a supercomputer. Each Iteration is 
made up of two stages and the red and black states succeed one another, 
however the black stage must not start until the preceding red stage has 
completed and conversely. 


W"~' 


51 


F*'». 




\3 


For the Dirichlet boundary conditions, if any term on the RHS 
in equation (6.3) and (6.4) belongs to the boundary, then the correspond- 
ing term is understood to have the prescribed boundary value. In case 
of the Newmann boundary condition, redefinition of the boundary data 
can be easily incorporated after each stage or two stages depending upon 
the number of grid points on the boundary, type of vector processor or 
tradeoff between the scalar and vector operations. 

It is worth noting that two-step Jacobi, which can use vector 
length equal to total number of nodes in the field, is also an attractive 
method for vector processor in which vector length is an important factor 
in the calculation rate and vector processor performance is at least 
twice its scalar performance. In many applications, the Jacobi method 
with acceleration parameter k = 1 may not be able to compete with the 
checkerboard SOR method. Frequently, cyclic change of red and black 
stages may give better convergence rate for the hopscotch method. 

There also exists a family of hopscotch methods such as line or 
zebra-like [44] and block methods [45, 46]. Some properties of hopscotch 
methods have been presented by Gourlay et. al. [47], Greenberg [48] 
employed the approximate factorization scheme of Beam and Warming [15] 
in hopscotch form and other hopscotch methods to investigate fluid 
dynamics problems on a serial computer. South et. al. [49] used Checker- 
board SOR, Zebra SOR and Checkerboard Leapfrog method for transonic 
flow calculations on the CDC-3TAR-100 and CRAY-1 vector computers. 

6 . 3 Implementation on the CYBER-203 

As the checkerboard SOR is the heart of an implicit method used in 
the present study, it was decided first to implement the model problem 
on the CYBER-203 and then to employ the same basic features 
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of the implementation for solving the Navier-Stokes equations. One 
possible way of approaching the model problem using the CYBER-200 
FORTRAN language will be discussed in this section. 

A close examination of the test problem indicates that the first 
task is to determine vectors of the red and black field variables from 
the arrays containing all field and boundary nodes of the same variables. 
Once the vectors of desired color are obtained, the arithmatic operations 
on these vectors are rather simple and can be performed using explicit 
vector instructions. An emphasis is made on the use of predefined 
vector functions and rich instruction set of the CYBER-200 FORTRAN 
compiler. The bit addressable memory, which allows the use of bit vectors 
is one of the important characterisitics of this machine. 

The total number of elements in any array equals to the product 
of its dimension sizes. All elements of the array are stored contiguously 
in a memory. To find the location of an array element for a given array 
T(A,B) of a particular instance of subscript T(a,c/, the formula 
a + A * (b-1) can be used. Thus an array can be thought of as a vector, 
and wherever required we will use word vector to represent an array in 
the remainder of this section. Each element of a bit vector requires 
storage of one bit in contrast to 64 bits required for each element of a 
single precision value vector (real or integer) on this machine. An 
element of a bit vector can take a val^e of either 1 or 0 representing 
logical operator truth and false respectively. All logical operations 
such as AND, NOT, etc. can be performed on the bit vectors. The logical 
operations are performed on either corresponding elements of two bit 
vector operands or a scalar operand paired with successive elements of 
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a vector operand. This important feature allows us to generate a bit 
vector of desired structure or pattern which in turn can be used as a 
control vector in some very efficient built-in functions. Also there 
are some functions which help to form an initial bit vector of some . 
desired 0-1 pattern. 

Some useful functions, which use bit vectors as control vectors will 
be briefly described, since they are an important part of the present 


implementation. Details of the builtin functions for the CYBER-200 
FORTRAN compiler can be found in reference [50], Bit vectors can be 
used as a control vector to select: elements from a value vector. The 
CMPRS function deletes selected elements from a real or integer vector 
as dictated by a bit control vector. The MERG function merges the 
elements in two value vectors into a result vector under control of a bit 
vector. The function CTRL changes the values of selected elements in a 
result vector using the values in an argument vector under the control 
of a bit control vector. 

For illustration, let us consider the model problem. As we would 
like to solve equation (6.3) in vector form, we must have vectors of 
all terms on the RHS of equation(6.3) . At the beginning of the first 
or red stage array f(5,5) of all nodes (including field and boundary) 
is available. The task is now to obtain vectors of all terms involving 
f and S on the RHS cf the equation. The procedure involves selecting 
and assembling all f^ n ’ r ^ at the red field nodes from vector f of the 
field arid boundary nodes. Similarly we must form vectors for 
f fj 11 ?^,, f[ n !^ , f , which are located at the black nodes. 

i+lj i» j~l i,j+l i, j-1 


Let us assume that all control bit vectors of desired pattern are stored 



in the memory and are available to facilitate the use of previously des- 
cribed functions. A 3 mentioned before, these bit vectors can be easily 
formed using some buit-in functions and logical operation such as 
AND, NOT, OR, etc. 

As shown in fig. 4 , execution of function CMPRS will give us the 
/ 

vector for f^ terms. Calling the CMPRS function with an appropriate 

control bit vector will give the vector for f^’^ (fig. 5). 

Similarly we can obtain vectors for the remaining f terms using appro- 

priate bit vectors. To obtain a vector for the source term S^ n ! r ^ , we 

can use the same bit vector as in fig. 4 , however the argument value 

vector will be S(5,5) instead of f(5,5). This completes the formation 

of all required vectors for solving equation (6.3) The equation involves 

scalar terms -j-, (1 - <) and h . These scalar terms are assumed to be 

implicitly expanded to the necessary vector length, with each element 

having the same scalar value. The arithmetic operations involved are 

straight-forward and equation (6.3) can be solved for f^ ’ on the 

*■9 J 

vector processor using explicit vector notations. 

Before we go to the second or black stage (eq. 6.4) we must update 
array f(5,5) using vector f^ n ^ ,r ^. One possible way to replace the 

+• 9 J 

old values of variable f at the red field nodes with the updated values, 
is to use the CTRL function. For the CTRL function, vector lengths of 
result, control bit and argument vector should be the same. The result 
vector in this case will be vector f and has the same vector length as 
of the control bit vector. However, the argument vector is about one 
half the length of vector f. Using the MERG function and a dummy 
value vector we can generate a vector of required length, having values 
of f. n !^ ,r ^ in the desired elements and the rest of the elements having 
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the values of the dummy vector. The values of the elements of a dummy 

vector are Insignificant and can be chosen arbitrarily. As shown in fig- 

6, MERG merges dummy vector DM, having arbitrary value for all elements 

and vector into a result vector RS, with the help of a control 

i»J 

vector. The merge stops when the result vector RS is full. Now having 
obtained vector RS of appropriate length and elements, the execution of 
the CTRL function with f as a result vector and RS as an argument vector 
under the control of a given bit vector will update the values of f 
at the red field nodes (fig. 7). This completes the implementation of 
the first or red stage (eq. 6.3) on the CYBER-203 computer using the bit 
control vector approach. 

The Implementation of the second or black stage (eq. 6. A) can be 
incorporated in the same framework using an updated array f and 
appropriate control bit vectors. Each iteration is made up of the two 
stages and the above iterative procedure for vector processing can be 
continued until desired accuracy is obtained. 

It is obvious that vector algorithms require more storage than 
scalar algorithms. However, due to large memory (1 million, 64 bit words) 
and sharing same storage locations the increased storage requirement can 
be handled properly in many applications. Instead of using the bit control 
vector approach, the above can be implemented using integer 

index vectors which are incorporated in functions such as GATHR (gather 
and SCATR (scatter). In these functions, instead of a bit vector, an 
integer vector with appropriate index values is used as a control vector. 
The integer index vector approach was not investigated in the present 
study. It Is also Interesting to note that in many applications the 
bit control vectors of desired structure need to be generated once only 
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and can be used many times in a computer code. Since each element of 
a bit vector corresponds to one bit only, the storage requirement for 
bit vectors are far less than conventional value vectors. 

6.4 The N a vier-Stokes Equations and Checkerboard SOR 

The governing equations for incompressible flow about an airfoil 
are the Navier-Stokes equations. In the present study, equations (3.21) 
and (3.22) for the velocities and equation (3.23) for the. pressure are 
solved simultaneously at each time step using the checkerboard SOR 
method. These transformed equations are somewhat complicated compared 
to the model problem and its implementation on the vector processor 
is more involved. 

The transformed or computational plane is rectangular regardless of 
the shape of the physical plane. The field nodes in the 2-D transformed 
plane are represented in a checkerboard pattern so that each red grid 
point has four black neighbors and vice versa. Three unknowns, two 
velocities and the pressure are associated with each node. All terms on 
the RHS of equations (3.21), (3.22) and (3.23) are represented using 
appropriate finite-difference approximations as discussed in Chapter V. 
Thus all ternm on the RHS of these difference equations can be represented 
in vector form as discussed for the test problem in the previous section. 
The storage of each term, including geometric coefficients, requires 
two vectors, each of about one-half the size of the entire field. 

Explicit vector instructions are employed to perform the arithmetic 
operations Involved in the equations. 

The transformed equations contain cross derivatives. The cross 
derivatives are evaluated using the central difference approximations. 

When solving the equations for a red field node, evaluation of cros3 
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derivatives involve red nodes in contradiction to the updated black 
nodes involved in evaluation of first and second derivatives. Although 
the cross derivatives are lagging by one stage per iteration in the 
present formulation, it did not show certain adverse effects on the 
convergence rate during numerical experimentation. Gourlay et. al. 

[51] has discussed handling of cross derivatives in some hopscotch 
methods . 

It is desirable to use the checkerboard SOR with relaxation parameter 
varying from iteration to iteration instead of a constant relaxation par- 
meter to accelerate convergence. The major unresolved problem concerning 
the checkerboard SOR is that of determination of sequence of optimum 
parameters which will produce the smallest number of iterations for a 
specified degree of convergence. For solution of the velocity equations 
using the classical SOR, the computation of sequence of acceleration 
parameter proposed by Thompson [52] and described in Appendix B pro- 
duced nea.ly optimal iterative procedure in previous investigations 
[29, 30]. In many cases the values of computed acceleration parameters 
for the checkerboard SOR and theoretical optimal acceleration parameters 
for the classical SOR are comparable and have about the same range 
[53]. The typical values of acceleration. parameters are less than one 
for the velocities and the pressure is not accelerated. The additional 
operations involved in computation of acceleration parameters is justified 
by producing faster convergence to the checkerboard SOR. 

The sequence used for solving the three governing equations simul- 
taneously may show some, if not impressive, improvement of the convergence 
rate. Out of several possible sequences, the sequence of solving 
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the velocity for the red nodes than the pressure for the black nodes 
again the velocity for the black nodes and the pressure for the red 
nodes is found to have favorable convergence characterisitcs. 

The details of all interesting features of the computer code and 
other studies will not be presented, partly due to the lack of space 
and partly because the outcome of some numerical experiments seems to 
be inconclusive. However, the present discussion shows that for the 
solution of large scientific problems on a vector computer, a consistent 
algorithm will always out perform an inconsistent algorithm implemented 
without considering the architecture of the computer. 
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Chapter VII 
COMPUTATIONAL RESULTS 

7 . 1 Coord i nate Systems 

Two different approaches discussed in Chapter II were used to 
generate "C" type coordinate systems for the NACA 66 ^- 01 ^ airofil, which 
is symmetric and has maximum thickness ratio of 18%. The grid contains 
II. points on the £ axis and JL points on the n axis, in particular, the 
values of II. and JL were set to 113 and 51 respectively for all coordi- 
nate systems. The major concern was to obtain accurate numerical reso- 
lution of the flow field about the airfoil. Since grid characteristics 
such as mesh spacing, smoothness and skewness can greatly affect the 
effectiveness of the hosted algorithm, it was decided to examine some 
effects of the grid characteristics on the flow field solution. 

It is desirable to have much finer grid spacing in the regions of 
boundary layers containing relatively high velocity gradient because a 
relatively coarse grid can lead to significant truncation errors in the 
solution of the Navier-Stokes equations. The RHS of the Poisson pressure 
equation contains velocity gradient terras and hence errors in dominant 
velocity gradient terms can result in erroneous values of the pressure 
near the body surface. The wall pressure boundary condition equation uses 
one sided difference approximation so errors in the pressure field near 
the wall can lead to errors in the implementation of the boundary condi- 
tions. The algebraic eddy viscosity turbulence model used in this study 
involves velocity gradient term and accurate computation of velocity 
gradients is important for consistent turbulent modeling. 
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Whenever a grid in the physical plane is not smooth the transforma- 
tion coefficients such as E , n , C and n can induce considerable 

X X’ y y 

numerical error in the solution caused by the nonunifora grid spacing. 

In some cases the grid skewness can also .lead to numerical oscillations 
and inaccuracies. A detailed discussion about the effects of these grid 
characteristics on the solution can be found in references [11,34], 

Three coordinate systems were used in the present effort. The first 
coordinate system C0RD1 (fig. 8) was generated using Sorenson's [35] 
approach and was rather crude. This coordinate system was used for 
development, testing and debugging of the computer code in the early 
stages. Coordinate systems C0RD2 (fig. 9) and C0RD3 (fig. 10) were 
generated using Thompson's [10] approach and Sorenson’s [35] approach 
respectively. The grid point distribution on the inner boundary was 
the same for these two coordinate systems. Also these two coordinate 
systems were extensively used for many numerical experiments and solu- 
tions to be presented in the remainder of this section. 

As central -difference approximations used in this study are suc- 
ceptible to numerical oscillations at higher Reynolds number, it was 
also decided to investiage effects of grid characteristics at lower 
Reynolds number to isolate oscillations caused by central-differences. 

Two coordinate systems COW32 and C0RD3, having the same grid point 
distribution on the inner boundary were tested for Reynolds number 1000. 
The solution is the trailing edge region had a dominant effect on the 
total flow field solution. Since the grid lines of C0RD2 (fig. 9) are 
skewed in the trailing edge region, coordinate systems C0RD3 (fig. 10) 
was generated with nearly orthogonal lines in the trailing edge region. 
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Better overall results were expected using coordinate system C0RD3, 
however it turned out the other way. Numerical results obtained using 
coordinate system C0RD2 were much better than those obtained using 
C0RD3. It was thought that other grid parameters such as coordinate 
stretching functions and the rate of change of grid spacing would have 
significant effect on the solution [34, 54}. An exponential stretching 
function was used in Sorensen's approach while in Thompson's approach co- 
ordinate control function was derived using the hyperbolic tangent as the 
point distribution function. As mentioned in reference [34], the hyper- 
bolic tangent is better than exponential and gives optimal truncation 
error. It should be noted here that the control function in Sorensen's 
approach controlled only spacing of the first line off the boundary 
and angle of inclination of £ * constant lines with the boundary while 
the control function in Thompson's approach was able to control grid 
line distribution (fig. 2) . The above case is not entirely conclusive, 
however it does show the importance of a proper coordinate system 
for the Navier-Stoke's solution. For a given problem finding of an 
optimum coordinate system by trial and error method is expensive, so 
we decided to limit our experimentation with coordinate systems. 

Some important parameters of two coordinate systems C0RD2 and 
C0RD3, which were used extensively in this study, are described below. 

For both grids the leading edge of the airfoil was located at (0,0) 
and the trailing edge was at (1,0). The y coordinates of the uppermost 
point at I 113 was +5.09 (5 chord lengths) the lowermost point at 
I = 1 was at -5.09. The x coordinate of the forward most point was 
-4.09 andof the backwardmost point was 11.0 (11 chord lengths). The 
value of index i at the leading edge was 57, The lower-surface trailing 
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edge and upper-surface trailing edge points were located at i = 21 and 
i * 93 respectively. The value of index i at the maximum airfoil thick- 
ness point on the lower and upper surface were 34 and 80 respectively , and 
the value of the x coordinate was 0.467. The term As denotes grid spac- 
ing between the first line off the boundary and boundary along £ =* con- 
stant lines. The values of As at the inner boundary for coordinate 
system C0RD2 were 0.000046, 0.000055, 0.000010 and 0.000026 at I = 1, 

21, 34 and 57 respectively. The minimum values of As at the inner 
boundary was 0.000001. The values of As at the. outer boundary were 
0.28, 0.35, 0.32 and 0.42 at I ■ 1, 21, 34, and 57 respectively. For 
the coordinate system C0RD3 the values of As at the inner boundary were 
set to 0.00001 and at the outer boundary were set to half the chord length 
(i.e. 0.5). The angles of inclination with which £ ■ constant li..es 
intersect the inner and outer boundaries were approximately 90°. For 
the two coordinate system the grid point distribution at the leading 
and trailing edge in the boundary layer region is shown in Table 1. 

Also note that coordinate system C0RD2 has a uniform spacing around the 
airfoil while C0RD3 has closer spacing at high curvature regions e.g. 
leading edge. 

7 . 2 Computational Procedure 

In this section a computational procedure for solving the govern- 
ing equations and some details of the comptuer code are described. . 

An appropriate coordinate system for a given problem was generated using 
separate grid generation computer code on a scalar machine. 

The values of x and y coordinates for a final grid were input to 
the. CYBER-200 FORTRAN compiler code. Bit control vectors of desired 
pattern were generated and stored. Since all geometric coefficients can 
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can be efficiently computed using bit vector approach, all tranformation 
coefficients were computed using bit vector approach, all transformation 
coefficients were computed by the computer code in the vector mode 
instead of using values of coefficients supplied by other scalar codes. 
Second-order accurate central-difference formulas were used to compute 
transformation coefficients in the field. On the airfoil surface upper 
down-stream boundary and lower down-stream boundary the transformation 
derviatves were computed using second order accurate one-sided forward 
or backward differences. All geometric coefficients were separated for 
the red and black nodes and stored. For restarting the flow from pre- 
viously obtained solutions, all required flow fiald variables and impor- 
tant parameters were read in. Before starting off a loop for time steps, 
all required bit control vectors for solution of the governing equations 
were generated and stored once and for all. 

All calculations were performed with the time step t = 0.01. The 
free-stream boundarycondition on theouter boundary was applied at every 
time step. The gradua start consisted of 100 time steps during which 
the free-stream velocities and the body force terms were given by 


u^ = t * for 0 < t < 1.0 (7.1) 

v m = t * g 2 (7.2) 

g-L - cos ip (7.3) 

g 2 = sin \p (7.4) 

after the gradual start 

u “ cos (7.5) 

°° for t > 1.0 


v = sin 

CO 

61 = 0 

6 2 a 0 


(7.6) 

(7.7) 

(7.8) 
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At every time step, two velocities and the pressure equations were 
solved simultaneously using checkerboard SOR. The sequence of solving 
the velocities for the red nodes, the pressure for the black nodes, the 
velocity for the black nodes and the pressure for the red nodes was used 
for each checkerboard iteration. The convergence criteria for each time 
step was established by the following procedure. The solution was either 
initially started or restarted from the previous time step. For the 
first time step first order two point backward-difference approximations 
were used for the time derivatives. The iteration continued until dif- 
ference between the magnitude of each flow variable (u, v and p) at two 
successive iterations were less than 0.0001. In most cases, maximum 
number of checkerboard iterations were limited to 30. At every iteration 
acceleration parameters for the velocities were computed using equations 
given in Appendix B. The computation of the accelerated parameters is 
more involved and requires considerable arithmetic operations. A flag 
was set when the solution converged within 10% of the established con- 
vergence criteria. When this flag was set, the computation of acclera- 
tion parameter was bypassed and iteration continued with the previously 
computed acceleration parameters to enhance the computational efficiency. 

The Neumann pressure boundary conditions, the re-entrant condition 

* 

and the downstream boundary condition was applied after every checker- 
board iteration. The trailing-edge pressure was extrapolated after 
applying the pressure boundary condition. This completes the solution 
procedure involved at every iteration. 

Once the solution converged within a given error norm or maximum 
number of iterations allowed were reached, the divergence of the velocity 
was computed .at every time step. As soon as the divergence of the 
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velocity was computed, it was smoothed out in most cases. Then the re- 
quired turbulence was switched on to compute eddy viscosity. Then the 
desired artificial viscosity was computed at every time step to incorpor- 
ate darning. A condition was established that artificial viscosity can- 
not be turned on unless turbulence was turned on. The above cycle was 
continued for the desired number of time steps. 

7 . 3 Some Numerical Experiments 

This section will present some numerical experiments carried out 
during the course of this study. It should be noted that all techniques 
described in this section were not tested thoroughly and some of them 
did not improve the solution or efficiency significantly. However many 
of the approaches attmpted, seem encouraging and may work well for other 
applications. The primary attention was focused upon the development of 
the efficient computer code and the computation of a reasonably accurate 
flow field solution using minimum computer resources. 

It is true, that the use of an appropriate algorithm is generally 
much more crucial than coding techniques. However, an optimized code , 
with proper algorithm can increase efficiency considerable. In this 
study once the algoritmh was settled upon, considerable time was spent on 
optimizing the code. There are very few loops in the code and they are 
generally unavoidable such as time step loop and iteration loop. Effort 
was concentrated to develop a code in the light of fundamental properties 
of the vector processor which allowed the use of explicit vector instruc- 
tions. All routines in the code were analyzed using a timing package 
which prints out a histogram of CPU usage. The main effort was diverted 
to some possible restructuring of routines and comparing its efficiency 
based on CPU timing. Initially all routines in the code could be compiled 
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with the highest level of optimization (B)) on the CYBER-200 FORTRAN com- 
piler. To incorporate various approaches to be tested, as discussed in 
the following paragraphs, a few routines were forced to one lower level 
of optimization (BE) . No attempts were made toward optimizing the memory 
and data raangaement procedures. Vectors of about 2720 length were employ- 
ed in the present study. Since the performance of the CYBER-203 in- 
creases with increase in vector length an application involving a very 
large coordinate system can result in relatively greater speed-up. Opti- 
mization may involve some work; however, for large scale problems usually 
it does pay off. 

Central differences used to approximate the spatial derivatives are 
easily succeptible to oscillations at higher Reynolds number. Computed 
solutions displayed large amplitude oscillations in the flow variables 
and destroyed the accuracy. The eddy viscosity model Increases the 
molecular viscosity and thus lowers the Reynolds number of the flow. The 
switching on of the turbulence model appeared to damp some oscillations 
but it did not show any significant degree of control over large ampli- 
tude oscillations. All attempts to obtain the steady state solution for 
flows at Reynolds number 10,000, which were started from fest, were 
unsuccesful beyond time t “ 1.0, even with inclusion of the turbulence 
model. One-sided differencing schemes may eliminate these nonlinear 
oscillations and may be vectorized from some simple regions. Hodge [25] 
used upwind differences for the first derivative terms, except the pres- 
sure gradient and velocity divergence terms, seemingly to avoid oscilla- 
tions caused by central-difference. For "C" and "0" type coordinate sys- 
tems, incorporation of the upwind-differences requrie checking the sign 

of contravariant velocity (E u + ( v) to incorporate appropriate indexing. 

x y 
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This condition may not allow efficient vectorization of one-sided differ- 
encing schemes for this study on the CYBER -203. 

One possible way to damp out the oscillations is by the use of an 
artificial viscosity. The adverse pressure gradient in the trailing edge 
regions had the dominant effect on the solution and it was assumed to 
trigger the nonlinear oscillations. The amplitude of these oscillations 
were small initially and remained localized near the trailing edge regions 
for some time, but in absence of damping its amplitude started increasing. 
The oscillations propagated toward the leading edge with passage of time. 
The flow was rather stable till time t ® 0.5, so in most cases it was de- 
cided to turn on damping at time t = 0.51, well before the oscillations 
started contaminating the solution. Several numerical experiments will be 
described before going into details of various forms of artifical viscosity. 

The time derivative of the divergence of velocity D^, appearing in 
equation (5.21) can be evaluated using either two point or three point 
backward difference approximations after the first time step. Several 
computer runs were made using both options. Error norms obtained using 
both cases for about the first one hundred time steps were almost the same 
indicating none was specifically responsible for divergence of the solu- 
tion. The spatial derivative terms of equation (5.21) were evaluated using 
second order accurate difference approxiamtions . Thus the use of two-point 
first order accurate approximation for the time derivaitve may reduce 
the overall accuracy of the equation. Since no emphasis was placed on the 
transient solution, and it was assumed that the time derivative term dis- 
appears in the steadt state, first order accurate approximation was used 
for the time derivative term of eq. (5.21) in most computer runs. Although 
it did not show any notic.able increase in computational efficiency, it is 
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interesting to note here that the two point backward approximations involve 
less computer operation and storage than the three point approximations. 


Since the implicit system of equations were solved at each time 
step by an iterative method, the previous time step solution was used 
as an initial guess for the next time step in all cases. For some cases 
a poor choice of initial guess may delay or destroy the convergence 
of the method. In an attempt to reduce the iterations by providing a 
good guess of the solution at the next time .level an initial guess which 
was close to the desired solution was tried. The initial guess for the 
velocities on the field and the re-entrant boundary was found using the 
following relations during the acceleration phase, 
n 


n-1 

u “ u + At cos 
ifj i» J 


v? j “ v 1 ) ! + At sin i|i 
i.j if j 


(7.9) 

(7.10) 


Instead of improving the convergence, the solution started diverg- 
ing. As the "C" type grid employed is coarse in the outer region and 
fine near the body and in the wake, the initial guess for the second 
attempt were found using the above relations only on the re-entrant 
section during the acceleration phase. Again no improvement 
was found and in both cases the solution started diverging approximately 
at time t = 0.5. For this study, it is not clear what should be the 
criteria to choose the initial guess in an effective manner and how 
it will accelerate the convergence. 

It was thought that the downstream boundary conditions may help 
control the oscillations or allov/ the passage of oscillations, which origi- 
nated in the trailing edge region. Instead of the downstream boundary 
conditions presented in section 4.2, the following downstream boundary 
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conditions were attempted. 
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The implementation of the free-stream boundary conditions, section 
A. 2, on the downstream boundary did not improve the solution during the 
acceleration phase. In another attempt, the velocity boundary condi- 
tions were the same as the free-stream boundary condition on the down- 
stream boundary, however the following pressure boundary condtion was 
used. 

P ^ = 0 (7 .14) 

Again, this boundary condition did not show any positive effect 
on the solution. Thus flow is perhaps much more sensitive to the outer 
and body surface boundary conditions with the downstream boundary con- 
dition having no significant influence on the solution. 

One possible way to enhance the stability of a numerical solution 
is to filter out unwanted oscillations using filters or smoothers. 

The use of smoothers will not eliminate the source of high amplitude 
oscillations but will control them by spreading them over some region. 
Since in some cases wavy solutions with high amplitude of flow quantites 
such as divergence of velocity, pressure can lead to unrealistic solution, 
the use of smoothers can help control oscillations. The divergence of 
velocity showed a wavy field in the direction o£ C = constant lines. 

Two types of divergence of velocity smoothers were attempted, one using 
the neighboring nodes in C a constant lines direction and the other 
using four neighboring nodes in both directions. Smoothing obtained 
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using four neighboring nodes was much better. It reduced the need of the 
artificial viscosity by some margin. The pressure field was wavy in 
the direction of n • constant lines. An attempt to smooth the pressure 
led to au incorrect solution. Next the source terms of the pressure 
equation (3.23) were smoothed oui . The solutions obtained using this 
approach were encouraging. Since the pressure equation was solved using 
an iterative method, for consistent smoothing the smoother should be 
applied at every iteration in contrast to the divergence smoother which 
was applied at every time step. The source term smoother may be compu- 
tationally inefficient due to computer operations and additional storage 
required. For some similar runs, results obtained using the source 
terra smoothers were about the same as those obtained using the divergenc 
smoother. Considering the above tests, the source term smoother was not 
practical way of smoothing the pressure oscillation in the present study, 
and hence the divergence smoother was employed for most computations. 

For turbulent flows, the values of eddy viscosity, computed using the two- 
layer algebraic model, varied considerably in the boundary layer and 
wake region. Some forms of artificial viscosities, to be presented 
later, were based on the eddy viscosity. Most of them used unsmoothed 
values of the eddy viscosity, however, a few of them employed smoothed 
values of the eddy viscosity. Artificial viscosities, which employed 
smoothed values of the eddy viscosity did not show any increase of 
effectiveness over the artificial viscosities based on unsmoothed eddy 
viscosity. . 

Perhaps the most effective way of eliminating the flow field 
oscillations caused by central-differencing at higher Reynolds number 



is to use an aritificial viscosity. In this study it was assumed that the 
flow was turbulent when an artificial viscosity was switched on and hence 
the molecular viscosity U *» 1 + e in the momentum equations was replaced 
by P =• 1 + £ + M a . Term U a denotes artificial viscosity and it increases 
the value of molecular viscosity. An artificial viscosity having uniform 
or constant value over the whole flow field will increase artificial 
diffusion everywhere in the field and is obviously not the solution of 
the problem. However, an artificial viscosity which is a function of some 
flow quantities having appreciable values in the region of extreme 
oscillations and negligible values everywhere else can effectively diffuse 
oscillations without changing characterisitcs of the original flow field 
considerably. Various forms of attempted artificial viscosities are 
listed below 
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Artificial viscosity was applied at every time step and results 
obtained using some of the above described forms of artificial viscosities 
will be presented in the next section. 
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7.4 Numerical Results 

The following general procedure was established for numerical compu- 
tation and it was common to many of the flow solutions attempted. A NACA 
66^018 airfoil section at zero angle of attack was considered for all 
computations. Coordinate systems with 113 (IL) grid points in £ direc- 
tion and 51 (JL) points in p direction were used. Gradual start was made 
up of 100 time steps with a time step size of 0.01. The previous time 
step solution was used as an initial guess for the next time step solu- 
tion. The acceleration parameters were computed using the local veloci- 

-4 

ties. The convergence criteria for the velocity and pressure were 10 
and the maximum number of iterations at each time step were limited to 
50. First order time differencing was used at the first time step and 

c 

second order time-differencing scheme was used for all subsequent time 
steps. The flow was laminar till time t = 0.5 . For turbulent flow, 
computation of eddy viscosity was turned on at time t = 0.51 and trans- 
ition was assumed to occur at minimum pressure on the upper and lower 
airfoil surfaces. For solutions with an artificial viscosity, the arti- 
ficial viscosity was turned on at time t “ 0.51 and it was computed at 
every time step. Also, whenever artificial viscosity was switched on, 
turbulence was assumed to be turned on at the same time. Except for some 
initial cases, the trailing edge pressure was extrapolated and the diver- 
gence of velocity smoother was used for the flow simulation. Some excep- 
tions to the above-described procedure will be mentioned at appropriate 
places in the following paragraphs. 

The first type of coordinate system C0RD1 considered in this study 
was generated using Sorenson's approach [35] and is shown in fig. 8. 

74 







The laminar flow past the airfoil was considered at Reynolds number 
of 10,000,. The pressure distributions at time t « 0.5 and t ■ 

0.7 are shown in fig. 11 and fig. 12 respectively. These pressure 
distributions indicate that the flow was well behaved until time t ® 

0., 7. The solution was restarted from time t =* 0.5 with the turbulence 
model switched on and transition occuring at the minimum pressure. 

The turbulent flow solution at time t = 0.7 was essentially the same 
as the laminar flow solution at the same time. At time t - 0.8 the 
pressure started oscillating at the trailing edge (fig. 13) . With the 
passage of time, the solution diverged at time t = 0.85. To damp out 
the trailing edge oscillations in the turbulent flow, the transition was 
forced to occur at the maximum airfoil thickness points on both surfaces 
instead of minimum pressure points which were almost at the trailing 
edge. As shown in fig. 14 the amplitude of the pressure oscillations 
was reduced somewhat at t ■ 0.8, however the solution diverged at 
time t * 0.89. Again turbulent flow was restarted from time t ■ 0.5 
with artificial viscosity added. The. artificial viscosity was computed 
using ftRej|V • Vj, with ft = 1.0. Previously observed oscillation 
disappeared at time t =* 0.8 [fig. 15] due to artificial diffusion. As 
expected, the pressure coefficient was going down with increase in time. 
Figure 16 shows the pressure distribution at time t =■ 1.0. With further 
increase in time, the solution diverged at t » 1.22. In pressure 
distributions for all stable solutions described so far, there was abrupt 
pressure rise at the trailing edge. To remedy this problem, the pressure 
at the trailing edge was extrapolated using eqs. (4.17-4.18). Fig. 17 shows 
the pressure distribuiton with the extrapolation for turbulent flow at 
t = 0.7. Another form of artificial viscosity, ftRej|v * V j , where 
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ft « e was attempted with transition occuring at maximum airfoil thick- 
ness points. The solution started oscillating at t = 1.0 (fig. 18). 

Once again an artificial viscosity , fiRejjV • v| with ft = 1.0, was 
attempted. However this time the artificial viscosity was computed at 
every iteration instead of every time step. Pressure distributions at 
time t = 1.0 and t = 2.0 are shown in fig. 19 and fig. 20 and the 
solution diverged at t = 2.28. Next an artificial viscosity 
ftRej/fw] [V • v| with ft = 1.0 was considered. Fig. 21 and fig 22 
show the pressure distribution at t ■=■ 1.7 and t « 2.0. At time t = 1.6 
the solution already started oscillating in the trailing edge region. 
Since coordinate system C0RD1 was rather crude, computations using 
C0RD1 were stopped. 

A second coordinate system C0RD2 (fig. 9) was generated using 
Thompson's approach [10] with control functions involving hyperbolic tan- 
gent to control n-line spacing in the boundary layer region. The Reynolds 
number considered was 10,000. The pressure distributions and leading - 
edge and trailing-edge velocity vectors at time t » 0.5 and t = 0.7 for 
the laminar flow are shown in fig. 23 and fig. 24. The turbulent flow 
was restarted from time t =* 0.5 with transition occuring at minimum 
pressure. Fig. 25 and fig. 26 shew the pressure distribution and 
velocity vectors at time t => 0.7 and t = 1.0. The laminar and turbu- 
lent solution at time t =* 0.7 were almost the same. An abrupt increase 
in the magnitude of velocity vectors at time t « 1.0 in the trailing 
edge region .indicates the presence of nonlinear oscillations in the 
solution. The turbulent flow solution diverged at t « 1.07. The 
turbulent flow with transistion forced at the maximum airfoil thick- 
ness was considered next and the solution at time t = 1.0 is shorn 
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in fig. 2.7. No significant improvement in the solution was found. 

An artificial viscosity, £2Rej|V • v| with fi «= e and transition occuring 
at maximum airfoil thickness was attempted. Comparing previous solutions 
with this solution at time t » 1.0 (fig. 28), no sufficient diffusion 
is the trailing edge region could be obtained. Perhaps, this was due 
to very small values of eddy viscosity which diluted the artificial 
viscosity. The solution diverged at time t = 1.32. An artificial 
viscosity fiReJ | V • v| with ft = 1.0 was considered. The turbulent flow 
solutions, using this artificial viscosity, at time t = 1.0 and t = 2.0 
are shown in fig. 29 and fig. 30. The solution diverged at time 
t = 2.18. These solutions indicated nonlinear oscillations in the 
trailing edge region, which could not be eliminated using the above 
described forms of artificial viscosity. It was thought that these 
oscillations were caused by skewed grid lines in the trailing edge 
region (fig. 9) which ultimately destroyed the solution. 

A third coordinate system C0RD3 (fig. 10) was generated with 
nearly orthogonal lines in the trailing edge region using Sorensen's 
approach. Laminar flow solutions for Reynolds number 10,000 at time 
t = 0.5 and t = 0.7 are shown in fig. 31 and fig. 32. The flow was 
well behaved as expected. Turbulent flow solution with transition 
occuring at minimum pressure was attempted and fig. 33 shows the 
solution at t = 0.9. Velocity vectors at the trailing edge reversed 
their direction with unusually large magnitude. This phenomena also 
indicates the presence of oscillations in the solution. The solution 
diverged at time t = 1.02. In an attempt to damp out these oscillations, 
transition was enforced at maximum airfoil thickness points. Fig. 34 
shows the pressure distribution and velocity vector plots at time 
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t: = 0.9. Note the magnitude of velocity vectors at the trailing edge 
has reduced somewhat but not sufficiently. At this stage, the divergence 
of velocity smoother and two-point backward differencing scheme for 
DD 

-“termof the pressure equation were incorporated. The solution was 
started from rest and the pressure distribution and velocity vectors 
at time t * 1.0 are shown in fig. 35. The RHS smoother, as discussed 
in Section 5.6 was also considered. The solution at time t = 1.0 is shown 
in fig. 36. No significant differences in the solutions obtained using 
these two smoothers were found. Smoothing of the pressure gave in- 
valid solutions. The smoothers were able to filter out some oscillations 
at the trailing edge. For all results presented hereafter, the diver- 
gence of velocity smoother was turned on and two-point backward dif- 

3D 

ferencing scheme for — term were employed all the time. Attention 

O U 

was now focused on the artificial viscosity. Fig. 37 shows the solution 
obtained at t = 1.0 using an artificial viscosity ftRejjv • v| with 

? 9 - 

3V 2 3V 2 ^ 

^ = {("TF") + (~V~) )*• Artificial viscosity with the same ft, but of 

ds of) 


form ftRej/faT] j"v • vj was considered next. Solution (Fig. 38) obtained 
using the latter form was somewhat better. Solutions obtained at t = 

1.0 using artificial viscosity ftRej|V • v| with ft = e|w| and ft = 1.0 
are shown in fig. 39 and fig. 40 respectively. An artificial viscosity 
based on Laplacian of the pressure i.e. ftReJAt|V P| with ft = 1.0 was 
attempted and the solution at t = 1.0 is shown in fig. 41. Artificial 
viscosity ftRej|V • v| with ft = [e^ e - 1.0] | to j and <j> = 4.0 gave some 
interesting results [fig. 42]. The values of 4> greater than 7.0 led 
to divergence of the solution for this particular case. The time 
history of solution obtained using artificial viscosity ftReJjv • v| 
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with & = <j>|w| and <j> = 1.0 is shown in figs. 43 - 46. The solution was 
almost steady at time t = 4.0 (fig. 46). This Reynolds number 10,000 
numerical solution was compared with Reynolds number 40,000 experimental 
solution [31] for qualitative purpose only. The discrepancies between 
the computational experiment results was thought due to unreasonably 
thick boundary layer and/or due to grid characteristics such as stretch- 
ing function. An attempt was made to restart flow from time t = 4.0 
without inclusion of the artificial viscosity, to obtain correct 
boundary layer. However, the solution diverged at time t = 4.28. 

Several values of <f>, ranging from 0.05 to 0.9 were experimented with but 
the results were not encouraging. Several computer runs with artificial 
viscosity given by ftRej|v • v| and ft = <J>e for values of <j> from 10 to 
1000 were made but without certain improvement. At this point, it was 
decided to lower the Reynolds number to isolate oscillations caused by 
higher Reynolds number and to investigate the effects of the coordinate 
systems on the solution. 

Coordinate system C0RD3 was used for a Reynolds number 1000. The 
pressure distribtuion and velocity vectors plots for the laminar flow 
at time t =* 1.0 and t = 1.5 are shown in fig. 47 and fig. 48. With 
increase in time the solution diverged at time 2.02 and the pressure 
distribution at time t = 2.0 is shown in fig. 49. Perhaps insufficient 
grid resolution in the boundary layer was responsible for the divergence 
of the solution. 

Next, coordinate system C0RD2. was considered for the laminar flow 
at a Reynold number of 1000. Time history of the flow at time t = 1.2, 
2.0, 3.0 and 4.0 is shorn in figs. 50 - 53. The flow characteristics 
were not changing with further increase in time and the number of 
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iterations came down to 8. Hence the flow was considered steady at 
time t= 4.0. For qualitative comparison, the numerical solution is 
compared with the experimental solution at a Reynolds number 40,000. 

The discrepancies in the pressure distribution from the leading edge 
to the maximum airfoil thickness can be identified. It is interesting 
to note here that coordinate system C0RD2 allowed to obtain the steady 
state solution while C0RD3 did not. A coordinate control function which 
was found using the hyperbolic tangent as the point distribution func- 
tion was used for coordinate system C0RD2. It was noted in reference 
[34] that a hyperbolic tangent function gave optimum truncation error. 

With the same coordinate system, i.e., C0RD2 an attempt was made to 
obtain Reynolds number 10,000 solution by restarting the laminar flow 
from Reynolds number 1J00 solution at t = 4.0. The pressure distribution 
and velocity vectors for Reynolds number 10,000 at t ** 5.0 is shown in 
fig. 54. 

It was thought that the accuracy of the solution during the 

acceleration phase had significant effect on the total flow field 

solution. Hence coordinate system C0RD2 was considered for Reynolds 

number 1000 laminar flow solution with increased number of iterations 

-4 

in initial stages. The error norm used was the same as before (10 ), 

however, the maximum number of iterations for initial time steps were 
increased to satisfy the above convergence criteria exactly. The 
number of iterations started increasing from 10 at the first time step 
to 75 at time t = 1.0, about 100 at time t = 1.5, about 90 at t = 2.0 
and about: 60 at t = 2.5. Some test runs were made with increased 
number of iterations and maximum number of ^iterations fixed at 50 beyond 
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time t = 3.0. No noticible difference between the solutions using 50 
and increased iterations was found. Hence the maximum number of iter- 
ations beyond time t = 3.0 were again fixed to 50. The time history 
of the solution starting with time t * 1.0 till time t = 10.0 is shown 
in figs. 55 - 64 at an interval of 100 time steps. Note the difference 
in the pressure distribution at t = 1.0 between this case (fig. 55) 
and a case with fixed iteration (fig. 50) . This difference in the 
pressure distribution becomes more obvious at time t = 2.0, 3.0 and 
4.0. For the present approach the solution had not achieved steady 
state at t = 4.0, with number of iteration about 40. Note that starting 
at time t = 3.0 velocity vectors at the leading edge region start 
changing its angle of inclination gradually and becoming parallel to 
the airfoil surfaces. Also the magnitude of velocity vectors in the 
trailing edge region keep increasing with passage of time. No notice- 
able difference in the pressure distribution was found between the 
solution at time t * 9.0 and t = 10.0 and henc.e computations for 
Reynolds number 1000 were stopped at t = 10.0. The pressure distribution 
in the leading edge region (fig. 64) has improved considerably compared 
to the previous steady state solution (fig. 53). Also, the experimental 
results at Reynolds number 40,000 matched qualitatively better than 
previous approaches. The pressure distribution in the leading edge 
region was the major cause of discrepancies. A closer look at coor- 
dinate system C0RD2 (fig. 9) shows that there is a sudden change in 
grid points spacing after approximately 9 points from the leading 
edge on both, upper and lower surfaces. Note that these points were 
placed by curvature of this surface. 
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Perhaps the redistribution of points in this region may help us to 
obtain correct pressure solution. 

Finally, Reynolds number was increased to 40,000 and the flow was 
restarted from the Reynolds number 1000 laminar solution at t = 10.0. 

The maximum number of iterations were limited to 50. The laminar flow 
solution diverged at t * 10.16 indicating the presence of large ampli- 
tude oscillations at higher Reynolds number. An attempt was made to 
damp out oscillations with the turbulence turned on at t » 10.01 and 
transition occuring at minimum pressure points. Again, the solution 
diverged at t » 10.17. Next the turbulent flow solution using artifi- 
cial viscosity ftRejjV • v| with ft “ <j>|w| and <j> <* 1.0 was considered. 

The use of the artificial viscosity allowed a steady solution to be 
obtained. Some minor oscillations In the trailing edge region were 
observed at about time t ■ 16.0. Hence the value of <}>was increased 
to 10.0 after time t ■ 16.5. The surface pressure distribution and 
the leading and trailing edge velocity vector plots at time t ** 20.0 
are shown in fig. 65. The separation was found to occur at about 60% 
chord position on the upper surface and at about 64.1% chord position 
on the lower surface. The computed surface pressure distribution Is 
compared with the experimental data. The use of artificial viscosity 
increased the thickness of the boundary layer. The surface grid point 
distribution was though to be responsible for the discrepancies between 
the computed and experimental surface pressure distribution in the lead- 
ing edge region. The distribution of the divergence of the velocity and 
the total viscosity (1 + e + p ) , for the first 20 n “ constant lines, at 
the leading and trailing edge is shown in Table 1. The maximum of the 
divergence of the velocity field occured at the second node point off 
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the wall at the leading edge and the value of the total viscosity at 
the same node was 1.11. Since the value of the eddy viscosity e was 
zero at the leading edge, the value of artificial viscosity was 0.11 at 
that point. With Increase in the value of J, the magnitude of the 
divergence started decreasing, however values of the aritifical viscos- 
ity was increasing till J = 10 and then it started decreasing. At the 
trailing edge the values of the divergence of the velocity were less 
compared to the values at the leading edge. However, the values of 
the artificial viscoisty at the trailing edge were larger than at the 
leading edge. The values were increasing with Increase in J till J » 

15 and then it started decreasing. The increase in the values of the 
artificial viscosity at the trailing edge was probably due to increase 
in the magnitude of vorticity and increase in the cell size. 

7 • 5 Computer tim e 

The computations were performed on the CDC, CYBER-203 computer at 
NASA Langley Research Center, Hampton, Virginia. For about: 50 iterations 
per time step, average CPU time for these computations was observed to 
be 3.3 seconds/ time step. This compares to 37.4 seconds/ time step [30] 
for a similar serial code on the CDC 7600. A factor of 11.36 was 
observed improvement in speed. A coordinate system with 5763 (113 x 
51) grid points was used in the present study, giving average computa- 

_5 

tional rate of 1.145 x 10 seconds/iteration/grid point. Further in- 
crease in speed through data management optimization and additional code 
seems possible. 
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Chapter VIII 
CONCLUSIONS 

The prime motivation of thi3 study was to develop a vectorizable 
algorithm for the implicit finite-difference solution of the incompres- 
sible Navier-Stokes equations in general curvlinear coordinates. The 
results indicate that it is economically feasible to obtain flow field 
soluiton past complex geometries. Much of the present effort was divert- 
ed to the numerical solution of the incompressible two-dimensional 
Reynolds averaged Navier-Stokes equations in nonconservative primitive 
variable formulation on the vector computer, especially to the development 
of a relaxation technique amenable to vector processing. The checkboard 
SOR relaxation technqiue and boundary-conforming coordinate system make 
the method efficient and versatile for a wide variety of configurations 
which could be addressed using a vector computer. The computer code was 
fully vectorized in the sense that all vectorizable loops were vectorized 
using explicit vector instructions and arithmetic operations were per- 
formed in a vector mode. The present computations on the CYBER-203 
indicated a speed gain of about 11 over CDC-7600. The acceleration para- 
meters, based on local velocities, were computed using the classical 
point SOR analysis and need to be studied in detail to make them optimal 
for the checkerboard SOR. The present implicit scheme is linearly un- 
conditionally stable excluding the pressure terms. This scheme with 
central-difference approximations for the spatial derivatives, exhibited 
oscillatory behavior at the higher Reynolds number. Out of several 
smoothers attempted, the divergence of velocity smoother proved to be an 
effective way of filtering oscillations during the early stages of the 
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solution. However, with passage of the time and increase amplitude of 
oscillations, the effectiveness of the smoother was lost and the accuracy 
of the solution was destroyed. Perhpas the use of an artificial viscosity 
was the most effective way of eliminating the flow field oscillations at 
higher Reynolds number for the present method. For solution at higher 
Reynolds number, restarting the flow from the steady state solution at 
lower Reynolds number was not particularly effective in controling the 
nonlinear oscillations. Thus initial conditions had little effect on the 
stability of the flow field solution. On the other hand, the accuracy 
of the solution during the accleration phase had significant influence 
on the steady solution. The down-stream boundary conditions showed little 
influence on the total flow field solution. Also it is not clear what 
type of initial guess for the checkerboard SOR can accelerate the conver- 
gence and reduce the number of iterations required for a given error norm. 
Computed results indicate that it is possible to obtain considerable speed- 
up using the present method. The effects of several coordinate systems 
on the numerical solution were studied. The importance of a proper coor- 
dinate line distribution to avoid grid induced errors and sensitiveness 
of the algorithm to the coordinate system were observed. In particular, 
the rate of change of line spacing in the boundary layer region was found 
to be more important than the grid line skewness at the boundary. Also, 
the grid point distribution on the body surface showed considerable in- 
fluence on the solution. The turbulence model used in this study shows 
little control over nonlinear oscillations. Hence some compromise in the 
flow field solution was made by using an artificial viscosity. The use 
of an artificial viscosity usually made the boundary layer thicker than 
what it should be, however it did stabilize the flow solution. The effects 
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of several forms of artificial viscosities on the solution were studied 
and compared for several test runs. Although steady state solutions 
were not attempted using each of them, the information about their 
relative merit can help to choose an appropriate form for particular 
application. Some ranges of a parameter, which was used to control 
the effect of various forms of airfoil viscosities, were obtained from 
numerical experiments. The use of various smoothers was found to reduce 
the need of an artificial viscosity by little margin. 

In the course of the present study the coordinate systems played 
a crucial role. The need for an optimized coordinate system for the 
Navier-Stokes solutions became apparent. Inability to compute flow 
field solution on one coordinate system and the discrepancies between the 
computed results and experimental data on the other coordinate system 
was perhaps due to the deficiencies of the coordinate systems. An 
adaptive coordinate system, which adapts to flow field variables 
gradients in a numerical solution may effectively solve this problem. 

The dynamic coupling of the coordinate governing equations with the 
flow field governing equations to resolve, developing gradients is 
perhaps the most promising approach to improve the overall outcome of 
the present computational procedure. An. adaptive grid may eliminate 
the extreme oscillations encountered using a fixed grid and enhance 
the performance of the algorithm because fewer iterations will be 
required due to improved convergence. Since computer operations become 
more efficient with increase in vector length on the CYBER-203 an 
application involving large number of grid points can increase the 
relative speed gain. The CYBER-203 is a one-pipe machine while the 
CYBER-205 is two -pipe machine. Thus significant speed up in computa- 
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tional rate can be obtained by running the computer code on the CYBER- 
205. The primitive variable formulation used for the governing equa- 
tions can be easily extended to three-dimensional problems. For large 
three-dimensional problems, the increased number of grid points saturate 
or nearly saturate the available memory. Hence some grid points mu3t 
be held in secondary storage and they must be transmitted to and from 
the. central memory. As overall execution time is a function of this 
memory transfer, the memory and data managment become much more important 
for three-dimensional problems. The computer time required for large 
scientific problem is generally so large that any increase in efficiency 
can represent substantial savings. In recent years, success in the 
development of high technology, such as very large scale integrated 
(VLSI) systems has revolutionized computer architecture. It seems 
possible to build a special purpose computer tailored for a special 
application. The software logic of a relaxation method, such as 
checkerboard SOR, can be implemented in hardware using VLSI elements. 

As point relaxation techniques are at the core of many systems of 
partial differential equations occuring in fluid dynamics, heat trans- 
fer, plasma dynamics, electrical network, semiconductor device model- 
ling and structural analysis, the hardware implementation cost can 
be justified. Success in the development of a special purpose computer 
in the future will significantly reduce real time processing and cost 
of computing. The future of high-speed scientific computing, with 
increased emphasis on vector processing, seems to be quite promising. 
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APPENDIX A 


VARIOUS RELATIONS AND DEFINITIONS IN THE TRANS FORMED PLANE 

This appendix contains the relations and definitions necessary to 
transform the equations of motion and boundary conditions from the physical 
to the computational plane. All transformations are presented in fully 
non-conservative form. The two following definitions are applicable 
throughout this appendix: 

f(x,y,t) = a scalar function with continuous first and second 
derivatives. 

F(x,y) * I F 1 (x,y) + j F 2 (x,y) = a vector function with con- 
tinuous first derivatives, i and j are Cartesian 
unit vectors. 

Definitions of the Transformation 


J ■ x _y - x.y 
n Z 3 n 

(A. 1) 

2 2 

a “ x + y 
n n 

(A. 2) 

8 « x.x + y.y 
C n Zn 

(A. 3) 

2,2 
Y “ x 5 + 

(A. 4) 

Dx - ax K - 28x^ + yx^ 

(A. 5) 

Dy “ cix_ r - 28x_ + yx 

ZZ nn 

(A. 6) 

0 “ (y^Dx - x^Dy) / J 

(A. 7) 

T ° (x Dy - y Dx) / J 
n n 

(A. 8) 
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Transformat ion of Scalar Derivatives 


f x - Of/9x) ytt - (y„f e - y 5 f n )-/J (a. 9) 

f y " < 3 f/ 3 y) x,t " ( Vn " W /J <a.io) 


f xx 

“ (3 2 f/3x 2 ) 

,t ■ (y n f EE 

• 2y t y n f tn + 

y E f nn )/j2 



+ (y 2 y - 

ky n y 55 

Vn y En + 

y E y nn ) ^ x n f £ ‘ 

' Vn )/j3 



. , 2 

+ ty„ X rr ~ 
n 55 

2y E y n*En + 

y S x m> ) (y E f n 

- 

(A. 11) 

f 

yy 

= (3 2 f/3y 2 )x 

■ c * <x 5 f « 

- 2x £ x n f En + 

x 2 f )J 2 

5 nn 



+ <v« - 

2x_x y_ + 

5 n'5n 

x 2 y ) (x f _ • 

5 nn n 5 

- V„ )/j3 



+ ( Vee - 

2x r x x, + : 
C n Cn 

£ 2 x ) (y_f ■ 

5 nn ' 3 E, n 

- y n f E )/j3 

(A. 12) 

f 

xy 

** (3 2 f/3x3y) 

t ° t( Vn + 

x y_)f_ - x 

n f 5 D i 

E y E f nn 



- V, f SE yji + [ V, X H - <x E y n + 

x n y E >x En 



+ x C y C x nn 1(y n f 5 ~ y S f n )/j3 + - l W« 

" ( Vn + Ve >y 6n + Wnn 1( Vn ’ x nV /j3 (A ' 13) 


f t - Of/at) x§y - Cf t ) ?>n - (f x :c t + f y y t ) (A. 14) 

Transformation of Vector Derivatives 
Laplacian: 

V 2 f - (af - 2Bf + yf ) / J 2 
55 5n nn 

- [ (Dx) (f ) + (Dy) (f )]/j 2 (A. 15) 

X y 


ORJGfNSIL r«,32 rs 

OF POOR quality 

V2f “ (af « * 23f ?n + vf nn + of n + Tf P /j2 


Gradient : 


- f " H V'. - Vn>! + ‘Vn ' V{ >jI/J 


Divergence; 


Normal to 5-line: 


? " ! 5/ I^J - (yi - x n j)/^r 

- fU T|-v 

Tangent to ri-line: 

t <n) = „<n) , . _ 

: “ x t - (x ? i + y? j)//7 


Tanget to 5-Lina: 


t (5) - n^> x k - - (x i + y j ) /a 
* »u ru 

ec ^ .tonal P er iva t ive s 

3f/3n (n) a n (n) • Vf » ( Y f n - Bf^/j/f 
3f/3t (r,) - t (n) . V f - f 7/f 

**’ V 

3f/3n <?) » n (?) . Vf - (of ? - By/j/^ 


(A. 16) 


(A. 17) 


' * - ’ " y 5 (F i ) n + W n “ x n (F 2 ^^ J (a. 18) 

■SSi-J * .Norma l. ..^id^-Tan^ent- Vectors 
Normal to n-line: 

n (n) - Vn/|Vn| - (-y £ i + x £ j)//7 (A>19) 


(A. 20) 


(A. 21) 


(A. 22) 


(A. 23) 


(A. 24) 


(A. 25) 





APPENDIX B 


The Navier -Stokes equations in non-conservative primitive variables 
formulation can be represented by the following general partial differen- 
tial equation, neglecting the cross derivative terms 


A,f-_ + A~f + B. f + B„f + Cf + D = 0 (B.l) 

i — nn i C ^ n 

where f denotes velocity u or v. For (IL-1) * (JL-1) simultaneous equations, 
the spectral radius 5^ of Jacobi iteration is given by [52]. 


|C - 2(A 1 + A 2 ) 


{ |4A^ - B^| co3(~^-) 


+ |4A^ - B^| cos(-j~-)} 


The optimal acceleration parameter k for SOR iterations can be obtained 


using the following relations 


1 + 1 - c' 


if 4A^ > and 4A^ > B^ (B.3) 


otherwise 


l + l + £ 


In this study, when conditions for equation (B.3) were satisfied, instead 
of computing acceleration parameter using equation (B.3) which gives 


k > 1.0, the acceleration parameters were set equal to 1. In other 
words, the momentum equations were solved using acceleration parameter 
less than or equal to 1.0. 






The coefficients in equation (B.l) are defined as follows. 




V = _E«_ 
1 „ ,2 
ReJ 


A„ - 
2 „ 2 
ReJ 


T 1 = ““2 + 2(x v -- y u) 
1 ReJ 2 n ^ 


T 2 = ~ 2(x_v - y u) 

ReJ 4 S 


For x momentum equation 


B i = T i + lb ( 2 V n “ Vn } 


b 2 = t 2 + isj ( - 2 v c + 


For y momentum equation 


= T i + uTr ^~ 2,j x + P y ) 

l 1 ReJ y n x-’n 
1 


B 2 = T 2 ReJ (2 ~ 

C = T 


At 


T = 1 for first-order time differencing 
3 

T = — for second-order time differencing 
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Coordinate System 
C0RD3 


Coordinate System 
C0RD2 


Coordinate System C0RD2 Re = 40,000 
Solution at t = 20.0 



I = 93 

I = 57 

1=93 

I = 57 

I = 

93 

I = 

57 


Arc length 

Arc length 

Arc length 

Arc length 

ABS (V • V) 

1 +£ + p 

ABS (V • V) 

1 + £ + p 

J 

S 

S 

S 

S 

' 

a 

■V — 


1 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2 

0.0000116 

0.0000116 

0.0000559 

0.0000026 

1.13 

3.75 

43.22 

1.07 

3 

0.0000269 

0.0000269 

0.0001273 

0.0000061 

1.54 

5.78 

50.96 

1.11 

4 

0.0000470 

0.0000470 

0.000220 

0.0000105 

1.53 

7.07 

50.41 

1.14 

5 

0.0000736 

0.0000734 

0.000341 

0.0000163 

1.15 

6.79 

50.25 

1.18 

6 

0.000108 

0.0000108 

0.000498 

0.0000239 

1.44 

10.18 

49.37 

1.23 

7 

0.000154 

0.000153 

0.000702 

0.0000336 

0.96 

8.76 

48.82 

1.29 

8 

0.000215 

0.000213 

0.000955 

0.000046 

1.25 

13.75 

47.48 

1.37 

9 

0.000295 

0.000292 

0.001309 

0.000062 

0.87 

12.04 

46.38 

1.46 

10 

0.00040 

0.000395 

0.001756 

0.000083 

1.03 

17.15 

44 .46 

1.57 

11 

0.000538 

0.000530 

0.002337 

0.000111 

0.80 

16.35 

42.73 

1.70 

12 

0.0007211 

0.000706 

0.003094 

0.000147 

0.78 

19.17 

40.15 

1.84 

13 

0.000962 

0.000935 

0.004079 ’ 

0.000193 

0.70 

20.44 

37.72 

2.01 

14 

0.001279 

0.00123 

0.005363 

0.000253 

0.54 

18.50 

34.46 

2.19 

15 

0.001697 

0.001621 

0.997039 

0.000331 

0.58 

22.48 

31.37 

2.38 

16 

0.00225 

0.002121 

0.00922 

0.000433 

0.33 

14.86 

27.53 

2.56 

17 

0.00297 

0.00276 

0.01206 

0.000565 

0.43 

20.52 

23.90 

2.73 

18 

0.00394 

0.00359 

0.01574 

0.000737 

0.10 

5.83 

19.67 

2.81 

19 

0.00521 

0.00464 

0.02049 

0.000962 

0.23 

13.22 

15.74 

2.84 

20 

0.00689 

0.00599 

0.02659 

0.00125 

0.11 

7.55 

11.39 

2.67 


Table 1. Coordinate Line Distribution and Divergence of Velocity and 
Artificial Viscosity Field at Leading and Trailing Edge. 
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OF POOR QUALITY 



Fig. 11 Surface Pressure Distribution 

Laminar Flow, t = 0.5, Re = 10,000 



Fig. 12 Surface Pressure. Distribution 

Laminar Flow, t = 0.7, Re = 10,000 
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Fig. 13 Surface Pressure Distribution 
Turbulent flow, t = 0.8, Re = 



Fig. 14 Surface Pressure Distribution 
Turbulent flow, t = 0,8, Re = 
Transition at Maximum Airfoil 
Points. 


10,000 


* 



* i * 

10,000 t. ’ 

Thickness * 

•! , 
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o .a A .i .a iJ3 

X 

Fig. 15 Surface Pressure Distribution 

Turbulent flow, t = 0.8, Re « 10,000 

V a ■=> QRe J 1 7 - V 1 , « =* 1.0 



Fig. 16 Surface. Pressure Distribution 

Turbulent flow, t = 1.0, Re » 10,000 

M “ ORe J I V • V I , $] *» 1.0 
a 1 „ , 1 




Fig. 17 Surface Pressure Distribution 

Turbulent flow, t = 0.7, Re » 10,000 
Extrapolated trailing edge pressure. 



Fig. 18 Surface Pressure Distribution 

Turbulent flow, t = 0.7, Re =■ 10,000 
Extrapolated trailing edge pressure 
M a ° fiReJ | V • V | , ft = e 



I i; 
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OF POOR QUALITY 



Fig. 19 Surface Pressure Distribution 

Laminar Flow, t = 1.0, Re = 10,000 
W a = f2Rej| V • V| , n ■ 1.0 

applied every iteration 



Fig. 20 Surface Pressure Distribution 

Laminar Flow, t = 2.0, Re = 10,000 
P a = fiReJ|v • V| , S2 « 1.0 

M a applied every iteration 
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Fig. 21 Surface Pressure Distribution 

Turbulent flow, t «• 1.6 Re ■=> .10,000 

u » fiRej/TIJT V • Vf, 0 ■ 1.0 
“ ^ « 



Fig. 22 Surface Pressure Distribution 

Turbulent f low, t = 2 .0, Re ■* 10,000 
U “ fiRej/Tur | V • Vl , h » 1.0 

3 v « « 
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(a) 



uc Eiu uinlui-uui Lit .n iiijil j n jiniihiniiiiil 
* .a .» .£ .t. ij) 

X 



(b) 



Fig. 23 


(a) Surface Pressure Distribution 

(b) Leading & Trailing Edge Velocity Vector Fields 
Laminar Flow, t ■ 0.5, Re = 10,000 
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Fig. 24 (a,' Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 0.7, Re * 10,000 
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(b) 


Fig. 25 


(a) Surface Pressure Distribution 

(b) Leading & Traillng-edge Velocity Vector Fields 
Turbulent Plow, t “ 0.7, Re » *10,000 
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(a) 





Fig. 26 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t » 1.0, Re = 10,000 
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Fig. 27 (a) Surface Pressure Distribution 

(b) Leading & Trail ing-edge Velocity Vector Fields 
Turbulent Flow, t = 1.0, Re = 10,000 
Transition at Maximum Airfoil Thickness Points 
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Fig. 28 (a) Surface Pressure Distribution | ■ , 

(b) Leading & Trailing-edge Velocity Vector Fields \ 

Turbulent Flow, t = 1.0, Re * 10,000 <> 

Transition at Maximum Airfoil Thickness Points 3 
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Fig. 29 (a) Surface Pressure Distribution 

(b) leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t = 1.0, Re = 10,000 

u a » nRej|v • v| , a = l.o 
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(a) 



(b) 



Fig. 30 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t » 2.0, Re = 10,000 
U = ORe J | V • V| , 0 = 1.0 
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Fig. 31 (a) Surface. Pressure Distributin 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 0.5, Re =» 10,000 
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(a) 






(a) Surface Pressure Distribution 

(b) Leading & Xrailing-edge Velocity Vector Fields 
Laminar flow, t =l 0.7, Re = 10,000 
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Fig. 33 (a) Surface Pressure Distribution 

(b) Leading & Trailing Edge Velocity Vector Fields 
Turbulent flow, t = 0.9, Re •> 10,000 
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Fig. 34 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent flow, t = 0.9, Re = 10,000 
Transition at Maximum Airfoil Thickness Points. 
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Fig. 35 (a) Surface Pressure Distribution j 

(b) Leading & Trailing-edge Velocity Vector Fields i 

Turbulent Flow, t = 1.0, Re “ 10,000 * 

Divergence of Velocity Smoothers, | 

1st Order Accurate 'i 

oC t 
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Fig. 36 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t = 1.0, Re ■» 10,000 
RHS Smoother 
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(b) 



Fig, 37 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t - 1.0, Re = 10,000 ^ 

p a = Rejjv • v| , « » {(8V 2 /35) 2 + (3V 2 /3n) 2 > 2 
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(b) 



fig. 38 


(a) 

Cb) 


Surface Pressure Distribution 
Leading & Trail ing-edge Velocity Vector 
Turbulent Fl ow, t ~ 1.0, Re = 10,000 
= SIReJ /) u] (V • VT, SI =» {(8V z /35) 2 + 


Fields 

(3V 2 /3n)¥* 
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Fig. 39 (a) Surface Pressure Distribution | 

(b) Leading & Trailing-edge Velocity Vector Fields !■ 

Turbulent Flow, t = 1.0, Re » 10,000 'J 

V, » LRejIV • Vi, ft = elwl ’,<■ 

® -V - r i 





Fig. 40 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t = 1.0, Re *» 10,000 
U a - fiRejjV • V| , 0 = 1.0 
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Fig. 41 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t “ 1.0, Re ■ 10,000 
U a « ftReJAt | V 2 P | , Q » 1.0 
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Fig. 42 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent flow, t “ 1.0, Re = 10,000 
Transition at Maximum Airfoil Thickness Points 
v = fiReJ | V • V| , SI = [| e - 1.0] | w l , <|> = 4.0 

A », «» 
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Fig. A3 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t = 1.0, Re “ 10,000 
H a ■ ftRe J | V • V|, S2»<J>|u>|, <j) » 1.0 
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Fig. 44 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t = 2.0, Re = 10,000 
u = fiRejIV • V j , ft = <j>|«j, (j> = 1.0 

a ^ «, ~ 
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Fig. 45 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent Flow, t = 3.0, Re =-- 10,000 
U a » ftRe J | V • V| , ii = <j>| w | , <j> = 1.0 
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Computational, Re «■ 10,000 

xxxx Experimental. Re - 40,000 





Fig., 46 (a) Surface Pressure Dtstribtuion 

(b) Leading & Trailing-edge Velocity Vector Fields 
Turbulent flow, t: = 4.0, Re = 10,000 
P a = ftReJ | V • V | , fi = '<> j u>| , <(>' » 1„0* 
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Fig. 47 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar flow, t a 1.0, Re = 1000 
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Fig. 48 ^a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Vector Velocity Plots 
Laminar Fbow, t =■ 1.5, Re “ 1000 
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Fig. 49 Surface Pressure Distribution 

Laminar Flow, t = 2.0, Re » 1000 
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Fig. 50 (a) Surface Pressure Distribtuion 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 1.0, Re “ 1000" 
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. Fig. 51 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 2.0, Re =* 1000. * 

t . 

. f 







142 







Fig. 52 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Veloctiy ’Vector Fields 
Laminar Flow, t. . = 3.0, Re ■-= 1000 
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— — Computational, Re ■ 1000 
xxxx Experimental, Re *> 40,000 
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Fig. 53 (a) Surface Pressure Distribution . 

(b) Leading & Trailing-edge Velocity. Vector Fields 
Laminar Flow, t >= 4.0, Re = 1000 
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Fig. 34 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 5.0, Re = 10,000 
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Fig. 57 (a) Surface Pressure Distribution < 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 3.0, Re “ 1000 
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Fig. 58 (a) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 4.0, Re = 1000 
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Fig. 60 


(a) Surface Pressure .Distr ibtuion - 

(b") Leading & Trailing-edge Velocity Vector Fields 
Laminar Floy, t “ 6.0, Re * lOpO 
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Fig. 61 (a) Surface Pressure Distribution 

(b)I.eading & irMling-edge VelociFy Vector Fields 
Laminar Flow, t = 7.0, Re “ 1000 . 
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Fig. 62 (a") Surface Pressure Distribution 

(b) Leading & Trailing -edge Velocity Vector Fields 
Laminar Flow, t: = 8.0, Re = 1000 • 
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Fig. 63 



(a) Surface Pressure Distribution 

(b) Leading 6 ^railing-edge Velocity Vector Fields 
Laminar Flow, t » 9.0, Re =.1000 . 
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Computational, Re ° 1000 

xxxx Experimental, Re " AO, 000 
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Fig. 6A (a*) Surface Pressure Distribution 

(b) Leading & Trailing-edge Velocity Vector Fields 
Laminar Flow, t = 10.0, Re =». 1000 . 
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(a) 



40,000 

40,000 


' (b) 



Fig. 65 (a) Surface Pressure Distribution, (b) Leading & Trailing 

Edge Velocity Vector Fields. Turbulent flow, t = 20.0, 
Re = 40,000. U a “fiRejjV • V J , <j)ju| 
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